Geant4  10.02.p02
G4LundStringFragmentation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4LundStringFragmentation.cc 91857 2015-08-07 13:55:49Z gcosmo $
28 // GEANT4 tag $Name: $ 1.8
29 //
30 // -----------------------------------------------------------------------------
31 // GEANT 4 class implementation file
32 //
33 // History: first implementation, Maxim Komogorov, 10-Jul-1998
34 // -----------------------------------------------------------------------------
35 #include "G4LundStringFragmentation.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "Randomize.hh"
39 #include "G4FragmentingString.hh"
40 #include "G4DiQuarks.hh"
41 #include "G4Quarks.hh"
42 
43 #include "G4Exp.hh"
44 #include "G4Pow.hh"
45 
46 //#define debug_LUNDfragmentation // Uzhi 20.06.2014
47 
48 // Class G4LundStringFragmentation
49 //*************************************************************************************
50 
51 G4LundStringFragmentation::G4LundStringFragmentation()
52 {
53 // ------ For estimation of a minimal string mass ---------------
54  Mass_of_light_quark =140.*MeV;
55  Mass_of_heavy_quark =500.*MeV;
56  Mass_of_string_junction=720.*MeV;
57 // ------ An estimated minimal string mass ----------------------
58  MinimalStringMass = 0.;
59  MinimalStringMass2 = 0.;
60 // ------ Minimal invariant mass used at a string fragmentation -
61  WminLUND = 0.45*GeV; //0.23*GeV; // Uzhi 0.7 -> 0.23 3.8.10 //0.8 1.5
62 // ------ Smooth parameter used at a string fragmentation for ---
63 // ------ smearing sharp mass cut-off ---------------------------
64  SmoothParam = 0.2;
65 
66  SetStringTensionParameter(1.);
67  SetDiquarkBreakProbability(0.05);
68  SetStrangenessSuppression(0.46); //(0.447); Uzhi 25.05.2015
69  SetDiquarkSuppression(0.05);
70 
71 // For treating of small string decays
72  for(G4int i=0; i<3; i++)
73  { for(G4int j=0; j<3; j++)
74  { for(G4int k=0; k<6; k++)
75  {
76  Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
77  }
78  }
79  }
80 //--------------------------
81  Meson[0][0][0]=111; // dbar-d Pi0
82  MesonWeight[0][0][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
83 
84  Meson[0][0][1]=221; // dbar-d Eta
85  MesonWeight[0][0][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
86 
87  Meson[0][0][2]=331; // dbar-d EtaPrime
88  MesonWeight[0][0][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
89 
90  Meson[0][0][3]=113; // dbar-d Rho0
91  MesonWeight[0][0][3]=pspin_meson*(1.-vectorMesonMix[0]);
92 
93  Meson[0][0][4]=223; // dbar-d Omega
94  MesonWeight[0][0][4]=pspin_meson*(vectorMesonMix[0]);
95 //--------------------------
96 
97  Meson[0][1][0]=211; // dbar-u Pi+
98  MesonWeight[0][1][0]=(1.-pspin_meson);
99 
100  Meson[0][1][1]=213; // dbar-u Rho+
101  MesonWeight[0][1][1]=pspin_meson;
102 //--------------------------
103 
104  Meson[0][2][0]=311; // dbar-s K0bar
105  MesonWeight[0][2][0]=(1.-pspin_meson);
106 
107  Meson[0][2][1]=313; // dbar-s K*0bar
108  MesonWeight[0][2][1]=pspin_meson;
109 //--------------------------
110 //--------------------------
111  Meson[1][0][0]=211; // ubar-d Pi-
112  MesonWeight[1][0][0]=(1.-pspin_meson);
113 
114  Meson[1][0][1]=213; // ubar-d Rho-
115  MesonWeight[1][0][1]=pspin_meson;
116 //--------------------------
117 
118  Meson[1][1][0]=111; // ubar-u Pi0
119  MesonWeight[1][1][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
120 
121  Meson[1][1][1]=221; // ubar-u Eta
122  MesonWeight[1][1][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
123 
124  Meson[1][1][2]=331; // ubar-u EtaPrime
125  MesonWeight[1][1][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
126 
127  Meson[1][1][3]=113; // ubar-u Rho0
128  MesonWeight[1][1][3]=pspin_meson*(1.-vectorMesonMix[0]);
129 
130  Meson[1][1][4]=223; // ubar-u Omega
131  //MesonWeight[1][1][4]=pspin_meson*(scalarMesonMix[0]);
132  MesonWeight[1][1][4]=pspin_meson*(vectorMesonMix[0]); // Uzhi 2015 scalar -> vector
133 //--------------------------
134 
135  Meson[1][2][0]=321; // ubar-s K-
136  MesonWeight[1][2][0]=(1.-pspin_meson);
137 
138  Meson[1][2][1]=323; // ubar-s K*-bar -
139  MesonWeight[1][2][1]=pspin_meson;
140 //--------------------------
141 //--------------------------
142 
143  Meson[2][0][0]=311; // sbar-d K0
144  MesonWeight[2][0][0]=(1.-pspin_meson);
145 
146  Meson[2][0][1]=313; // sbar-d K*0
147  MesonWeight[2][0][1]=pspin_meson;
148 //--------------------------
149 
150  Meson[2][1][0]=321; // sbar-u K+
151  MesonWeight[2][1][0]=(1.-pspin_meson);
152 
153  Meson[2][1][1]=323; // sbar-u K*+
154  MesonWeight[2][1][1]=pspin_meson;
155 //--------------------------
156 
157  Meson[2][2][0]=221; // sbar-s Eta
158  MesonWeight[2][2][0]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
159 
160  Meson[2][2][1]=331; // sbar-s EtaPrime
161  MesonWeight[2][2][1]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
162 
163  Meson[2][2][3]=333; // sbar-s EtaPrime
164  MesonWeight[2][2][3]=pspin_meson*(vectorMesonMix[5]);
165 //--------------------------
166 
167  for(G4int i=0; i<3; i++)
168  { for(G4int j=0; j<3; j++)
169  { for(G4int k=0; k<3; k++)
170  { for(G4int l=0; l<4; l++)
171  { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
172  }
173  }
174  }
175 
176  G4double pspin_barion_in=pspin_barion;
177  //pspin_barion=0.75;
178 //---------------------------------------
179  Baryon[0][0][0][0]=1114; // Delta-
180  BaryonWeight[0][0][0][0]=1.;
181 
182 //---------------------------------------
183  Baryon[0][0][1][0]=2112; // neutron
184  BaryonWeight[0][0][1][0]=1.-pspin_barion;
185 
186  Baryon[0][0][1][1]=2114; // Delta0
187  BaryonWeight[0][0][1][1]=pspin_barion;
188 
189 //---------------------------------------
190  Baryon[0][0][2][0]=3112; // Sigma-
191  BaryonWeight[0][0][2][0]=1.-pspin_barion;
192 
193  Baryon[0][0][2][1]=3114; // Sigma*-
194  BaryonWeight[0][0][2][1]=pspin_barion;
195 
196 //---------------------------------------
197  Baryon[0][1][0][0]=2112; // neutron
198  BaryonWeight[0][1][0][0]=1.-pspin_barion;
199 
200  Baryon[0][1][0][1]=2114; // Delta0
201  BaryonWeight[0][1][0][1]=pspin_barion;
202 
203 //---------------------------------------
204  Baryon[0][1][1][0]=2212; // proton
205  BaryonWeight[0][1][1][0]=1.-pspin_barion;
206 
207  Baryon[0][1][1][1]=2214; // Delta+
208  BaryonWeight[0][1][1][1]=pspin_barion;
209 
210 //---------------------------------------
211  Baryon[0][1][2][0]=3122; // Lambda
212  BaryonWeight[0][1][2][0]=(1.-pspin_barion)*0.5;
213 
214  Baryon[0][1][2][1]=3212; // Sigma0
215  BaryonWeight[0][1][2][1]=(1.-pspin_barion)*0.5;
216 
217  Baryon[0][1][2][2]=3214; // Sigma*0
218  BaryonWeight[0][1][2][2]=pspin_barion;
219 
220 //---------------------------------------
221  Baryon[0][2][0][0]=3112; // Sigma-
222  BaryonWeight[0][2][0][0]=1.-pspin_barion;
223 
224  Baryon[0][2][0][1]=3114; // Sigma*-
225  BaryonWeight[0][2][0][1]=pspin_barion;
226 
227 //---------------------------------------
228  Baryon[0][2][1][0]=3122; // Lambda
229  BaryonWeight[0][2][1][0]=(1.-pspin_barion)*0.5;
230 
231  Baryon[0][2][1][1]=3212; // Sigma0
232  BaryonWeight[0][2][1][1]=(1.-pspin_barion)*0.5;
233 
234  Baryon[0][2][1][2]=3214; // Sigma*0
235  BaryonWeight[0][2][1][2]=pspin_barion;
236 
237 //---------------------------------------
238  Baryon[0][2][2][0]=3312; // Theta-
239  BaryonWeight[0][2][2][0]=1.-pspin_barion;
240 
241  Baryon[0][2][2][1]=3314; // Theta*-
242  BaryonWeight[0][2][2][1]=pspin_barion;
243 
244 //---------------------------------------
245 //---------------------------------------
246  Baryon[1][0][0][0]=2112; // neutron
247  BaryonWeight[1][0][0][0]=1.-pspin_barion;
248 
249  Baryon[1][0][0][1]=2114; // Delta0
250  BaryonWeight[1][0][0][1]=pspin_barion;
251 
252 //---------------------------------------
253  Baryon[1][0][1][0]=2212; // proton
254  BaryonWeight[1][0][1][0]=1.-pspin_barion;
255 
256  Baryon[1][0][1][1]=2214; // Delta+
257  BaryonWeight[1][0][1][1]=pspin_barion;
258 
259 //---------------------------------------
260  Baryon[1][0][2][0]=3122; // Lambda
261  BaryonWeight[1][0][2][0]=(1.-pspin_barion)*0.5;
262 
263  Baryon[1][0][2][1]=3212; // Sigma0
264  BaryonWeight[1][0][2][1]=(1.-pspin_barion)*0.5;
265 
266  Baryon[1][0][2][2]=3214; // Sigma*0
267  BaryonWeight[1][0][2][2]=pspin_barion;
268 
269 //---------------------------------------
270  Baryon[1][1][0][0]=2212; // proton
271  BaryonWeight[1][1][0][0]=1.-pspin_barion;
272 
273  Baryon[1][1][0][1]=2214; // Delta+
274  BaryonWeight[1][1][0][1]=pspin_barion;
275 
276 //---------------------------------------
277  Baryon[1][1][1][0]=2224; // Delta++
278  BaryonWeight[1][1][1][0]=1.;
279 
280 //---------------------------------------
281  Baryon[1][1][2][0]=3222; // Sigma+
282  BaryonWeight[1][1][2][0]=1.-pspin_barion;
283 
284  Baryon[1][1][2][1]=3224; // Sigma*+
285  BaryonWeight[1][1][2][1]=pspin_barion;
286 
287 //---------------------------------------
288  Baryon[1][2][0][0]=3122; // Lambda
289  BaryonWeight[1][2][0][0]=(1.-pspin_barion)*0.5;
290 
291  Baryon[1][2][0][1]=3212; // Sigma0
292  BaryonWeight[1][2][0][1]=(1.-pspin_barion)*0.5;
293 
294  Baryon[1][2][0][2]=3214; // Sigma*0
295  BaryonWeight[1][2][0][2]=pspin_barion;
296 
297 //---------------------------------------
298  Baryon[1][2][1][0]=3222; // Sigma+
299  BaryonWeight[1][2][1][0]=1.-pspin_barion;
300 
301  Baryon[1][2][1][1]=3224; // Sigma*+
302  BaryonWeight[1][2][1][1]=pspin_barion;
303 
304 //---------------------------------------
305  Baryon[1][2][2][0]=3322; // Theta0
306  BaryonWeight[1][2][2][0]=1.-pspin_barion;
307 
308  Baryon[1][2][2][1]=3324; // Theta*0
309  BaryonWeight[1][2][2][1]=pspin_barion;
310 
311 //---------------------------------------
312 //---------------------------------------
313  Baryon[2][0][0][0]=3112; // Sigma-
314  BaryonWeight[2][0][0][0]=1.-pspin_barion;
315 
316  Baryon[2][0][0][1]=3114; // Sigma*-
317  BaryonWeight[2][0][0][1]=pspin_barion;
318 
319 //---------------------------------------
320  Baryon[2][0][1][0]=3122; // Lambda
321  BaryonWeight[2][0][1][0]=(1.-pspin_barion)*0.5;
322 
323  Baryon[2][0][1][1]=3212; // Sigma0
324  BaryonWeight[2][0][1][1]=(1.-pspin_barion)*0.5;
325 
326  Baryon[2][0][1][2]=3214; // Sigma*0
327  BaryonWeight[2][0][1][2]=pspin_barion;
328 
329 //---------------------------------------
330  Baryon[2][0][2][0]=3312; // Sigma-
331  BaryonWeight[2][0][2][0]=1.-pspin_barion;
332 
333  Baryon[2][0][2][1]=3314; // Sigma*-
334  BaryonWeight[2][0][2][1]=pspin_barion;
335 
336 //---------------------------------------
337  Baryon[2][1][0][0]=3122; // Lambda
338  BaryonWeight[2][1][0][0]=(1.-pspin_barion)*0.5;
339 
340  Baryon[2][1][0][1]=3212; // Sigma0
341  BaryonWeight[2][1][0][1]=(1.-pspin_barion)*0.5;
342 
343  Baryon[2][1][0][2]=3214; // Sigma*0
344  BaryonWeight[2][1][0][2]=pspin_barion;
345 
346 //---------------------------------------
347  Baryon[2][1][1][0]=3222; // Sigma+
348  BaryonWeight[2][1][1][0]=1.-pspin_barion;
349 
350  Baryon[2][1][1][1]=3224; // Sigma*+
351  BaryonWeight[2][1][1][1]=pspin_barion;
352 
353 //---------------------------------------
354  Baryon[2][1][2][0]=3322; // Theta0
355  BaryonWeight[2][1][2][0]=1.-pspin_barion;
356 
357  Baryon[2][1][2][1]=3324; // Theta*0
358  BaryonWeight[2][1][2][1]=pspin_barion;
359 
360 //---------------------------------------
361  Baryon[2][2][0][0]=3312; // Theta-
362  BaryonWeight[2][2][0][0]=1.-pspin_barion;
363 
364  Baryon[2][2][0][1]=3314; // Theta*-
365  BaryonWeight[2][2][0][1]=pspin_barion;
366 
367 //---------------------------------------
368  Baryon[2][2][1][0]=3322; // Theta0
369  BaryonWeight[2][2][1][0]=1.-pspin_barion;
370 
371  Baryon[2][2][1][1]=3324; // Theta*0
372  BaryonWeight[2][2][1][1]=pspin_barion;
373 
374 //---------------------------------------
375  Baryon[2][2][2][0]=3334; // Omega
376  BaryonWeight[2][2][2][0]=1.;
377 
378 //---------------------------------------
379  pspin_barion=pspin_barion_in;
380  /*
381  for(G4int i=0; i<3; i++)
382  { for(G4int j=0; j<3; j++)
383  { for(G4int k=0; k<3; k++)
384  { for(G4int l=0; l<4; l++)
385  { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;}
386  }
387  }
388  }
389  G4int Uzhi;
390  G4cin>>Uzhi;
391  */
392 
393  SetStrangenessSuppression(0.375); // Uzhi May 2015
394  Prob_QQbar[0]=StrangeSuppress; // Probability of ddbar production
395  Prob_QQbar[1]=StrangeSuppress; // Probability of uubar production
396  Prob_QQbar[2]=1.0-2.*StrangeSuppress; // Probability of ssbar production
397  SetStrangenessSuppression(0.46); //(0.447); // Uzhi May 2014
398 
399  //A.R. 25-Jul-2012 : Coverity fix.
400  for ( G4int i=0 ; i<35 ; i++ ) {
401  FS_LeftHadron[i] = 0;
402  FS_RightHadron[i] = 0;
403  FS_Weight[i] = 0.0;
404  }
405  NumberOf_FS = 0;
406 
407 }
408 
409 // --------------------------------------------------------------
410 G4LundStringFragmentation::~G4LundStringFragmentation()
411 {}
412 
413 
414 //--------------------------------------------------------------------------------------
415 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
416 {
417  G4double EstimatedMass=0.;
418  G4int Number_of_quarks=0;
419  G4int Number_of_squarks=0;
420 
421  G4double StringM=string->Get4Momentum().mag();
422 
423  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
424 
425 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
426 // G4cout<<"MinStringMass// Input String mass "<<string->Get4Momentum().mag()<<" Qleft "<<Qleft<<G4endl;
427 #endif
428 
429  if( Qleft > 1000)
430  {
431  Number_of_quarks+=2;
432  G4int q1=Qleft/1000;
433  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
434  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
435 
436  G4int q2=(Qleft/100)%10;
437  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
438  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
439 // EstimatedMass +=Mass_of_string_junction;
440 //if((q1 > 2)||(q2 > 2)) EstimatedMass -= 120*MeV; // Uzhi April 2014 ???
441  }
442  else
443  {
444  Number_of_quarks++;
445  if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
446  if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
447  }
448 
449 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
450 // G4cout<<"Min mass with Qleft "<<Qleft<<" "<<EstimatedMass<<G4endl;
451 #endif
452  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
453  if( Qright > 1000)
454  {
455  Number_of_quarks+=2;
456  G4int q1=Qright/1000;
457  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
458  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
459 
460  G4int q2=(Qright/100)%10;
461  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
462  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
463  //EstimatedMass +=Mass_of_string_junction;
464 //if((q1 > 2)||(q2 > 2)) EstimatedMass -= 120*MeV; // Uzhi April 2014 ???
465  }
466  else
467  {
468  Number_of_quarks++;
469  if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
470  if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
471  }
472 
473 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
474 // G4cout<<"Min mass with Qleft and Qright "<<Qright<<" "<<EstimatedMass<<G4endl;
475 // G4cout<<"Number_of_quarks "<<Number_of_quarks<<" Number_of_squarks "<<Number_of_squarks<<G4endl;
476 #endif
477 
478  if(Number_of_quarks==2){EstimatedMass += 70.*MeV;} //100.*MeV;}
479 // if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
480  if(Number_of_quarks==3)
481  {
482  if(Number_of_squarks==0) {EstimatedMass += 740.*MeV;} // 700 Uzhi July 2014
483  if(Number_of_squarks==1) {EstimatedMass += 740.*MeV;} // 740 Uzhi Nov 2014
484  if(Number_of_squarks==2) {EstimatedMass += 400.*MeV;}
485  if(Number_of_squarks==3) {EstimatedMass += 382.*MeV;}
486  }
487  if(Number_of_quarks==4)
488  {
489  if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2020.;}//1880.;}
490  // if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2051.;}
491  else if((StringM > 2232.) && ( EstimatedMass < 2730)){EstimatedMass = 2570.;}
492  else if((StringM > 5130.) && ( EstimatedMass < 3450)){EstimatedMass = 5130.;}
493  else
494  {
495 // VU 30 May 2014 EstimatedMass -=2.*Mass_of_string_junction;
496  if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
497  else {EstimatedMass+=100.*MeV;}
498  }
499  }
500 
501 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
502 // G4cout<<"EstimatedMass -------------------- "<<EstimatedMass <<G4endl;
503 #endif
504  MinimalStringMass=EstimatedMass;
505  SetMinimalStringMass2(EstimatedMass);
506 }
507 
508 //--------------------------------------------------------------------------------------
509 void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue)
510 {
511  MinimalStringMass2=aValue * aValue;
512 }
513 
514 //--------------------------------------------------------------------------------------
515 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
516 {
517  // Can no longer modify Parameters for Fragmentation.
518  PastInitPhase=true;
519 
520  SetMassCut(160.*MeV); // For LightFragmentationTest it is required
521  // that no one pi-meson can be produced.
522 
523  G4FragmentingString aString(theString);
524  SetMinimalStringMass(&aString);
525 
526 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
527  G4cout<<G4endl<<"LUND StringFragm: String Mass "
528  <<theString.Get4Momentum().mag()<<" Pz "
529  <<theString.Get4Momentum().pz()
530  <<"------------------------------------"<<G4endl;
531  G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
532  <<theString.GetRightParton()->GetPDGcode()<<" "
533  <<theString.GetDirection()<< G4endl;
534  G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
535  G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
536  G4cout<<"Check for Fragmentation "<<G4endl;
537 #endif
538 
539  G4KineticTrackVector * LeftVector(0);
540 
541 //Uzhi 29.05.2015 if(!IsFragmentable(&aString)) // produce 1 hadron True ===============
542  if(!aString.FourQuarkString() && !IsFragmentable(&aString))
543  {
544 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
545  G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
546 #endif
547  SetMassCut(10000.*MeV);
548  LeftVector=LightFragmentationTest(&theString);
549  SetMassCut(160.*MeV);
550 
551  LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
552  LeftVector->operator[](0)->SetPosition(theString.GetPosition());
553 
554  if(LeftVector->size() > 1)
555  {
556  // 2 hadrons created from qq-qqbar are stored
557  LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
558  LeftVector->operator[](1)->SetPosition(theString.GetPosition());
559  }
560  return LeftVector;
561  } // end of if(!IsFragmentable(&aString))
562 
563 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
564  G4cout<<"The string will be fragmented. "<<G4endl;
565 #endif
566 
567  // The string can fragment. At least two particles can be produced.
568  LeftVector =new G4KineticTrackVector;
569  G4KineticTrackVector * RightVector=new G4KineticTrackVector;
570 
571  G4ExcitedString *theStringInCMS=CPExcited(theString);
572 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
573  G4cout<<"CMS Left mom "<<theStringInCMS->GetLeftParton()->Get4Momentum()<<G4endl;
574  G4cout<<"CMS Right mom "<<theStringInCMS->GetRightParton()->Get4Momentum()<<G4endl;
575 #endif
576  G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
577 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
578  G4cout<<"aligCMS Left mom "<<theStringInCMS->GetLeftParton()->Get4Momentum()<<G4endl;
579  G4cout<<"aligCMS Right mom "<<theStringInCMS->GetRightParton()->Get4Momentum()<<G4endl;
580  G4cout<<G4endl<<"LUND StringFragm: String Mass "
581  <<theStringInCMS->Get4Momentum().mag()<<" Pz Lab "
582  <<theStringInCMS->Get4Momentum().pz()
583  <<"------------------------------------"<<G4endl;
584  G4cout<<"String ends and Direction "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
585  <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
586  <<theStringInCMS->GetDirection()<< G4endl;
587 #endif
588  G4bool success = Loop_toFragmentString(theStringInCMS, LeftVector, RightVector);
589 
590  delete theStringInCMS;
591 
592  if ( ! success )
593  {
594  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
595  LeftVector->clear();
596  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
597  delete RightVector;
598  return LeftVector;
599  }
600 
601  // Join Left- and RightVector into LeftVector in correct order.
602  while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
603  {
604  LeftVector->push_back(RightVector->back());
605  RightVector->erase(RightVector->end()-1);
606  }
607  delete RightVector;
608 
609  CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
610 
611  G4LorentzRotation toObserverFrame(toCms.inverse());
612 
613  G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
614  G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
615 
616  for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
617  {
618  G4KineticTrack* Hadron = LeftVector->operator[](C1);
619  G4LorentzVector Momentum = Hadron->Get4Momentum();
620  //G4cout<<"Hadron "<<Hadron->GetDefinition()->GetParticleName()<<" "<<Momentum<<G4endl;
621  Momentum = toObserverFrame*Momentum;
622  Hadron->Set4Momentum(Momentum);
623 
624  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
625  Momentum = toObserverFrame*Coordinate;
626  Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
627  G4ThreeVector aPosition(Momentum.vect());
628  Hadron->SetPosition(PositionOftheStringCreation+aPosition);
629  };
630 
631  return LeftVector;
632 }
633 
634 //----------------------------------------------------------------------------------
635 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
636 {
637  SetMinimalStringMass(string);
638 // return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
639 // G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
640  return MinimalStringMass < string->Get4Momentum().mag();
641 }
642 
643 //----------------------------------------------------------------------------------------
644 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
645 {
646  SetMinimalStringMass(string);
647 
648  if (string->FourQuarkString())
649  {
650  return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
651  } else {
652  // Uzhi 23 Jan. 2015 0.88 -> 0.66
653  return G4UniformRand() < G4Exp(-0.66e-6*(string->Mass()*string->Mass() -
654  MinimalStringMass*MinimalStringMass));
655  }
656 }
657 
658 //----------------------------------------------------------------------------------------------------------
659 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
660  G4KineticTrackVector * LeftVector,
661  G4KineticTrackVector * RightVector)
662 {
663  //... perform last cluster decay
664 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
665  G4cout<<"Split last-----------------------------------------"<<G4endl;
666 #endif
667  G4LorentzVector Str4Mom=string->Get4Momentum();
668  G4ThreeVector ClusterVel=string->Get4Momentum().boostVector();
669  G4double StringMass=string->Mass();
670 
671  G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
672 
673  NumberOf_FS=0;
674  for(G4int i=0; i<35; i++) {FS_Weight[i]=0.;}
675 
676 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
677  G4cout<<"StrMass "<<StringMass<<" q "
678  <<string->GetLeftParton()->GetParticleName()<<" "
679  <<string->GetRightParton()->GetParticleName()<<G4endl;
680 #endif
681 
682  string->SetLeftPartonStable(); // to query quark contents..
683 
684  if (string->FourQuarkString() )
685  {
686  // The string is qq-qqbar type. Diquarks are on the string ends
687  //G4cout<<"The string is qq-qqbar type. Diquarks are on the string ends"<<G4endl;
688 
689  if(StringMass-MinimalStringMass < 0.)
690  {
691  if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
692  {
693  return false;
694  }
695  } else
696  {
697  Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
698 
699  if(NumberOf_FS == 0) return false;
700 
701  G4int sampledState = SampleState();
702  if(string->GetLeftParton()->GetPDGEncoding() < 0)
703  {
704  LeftHadron =FS_LeftHadron[sampledState];
705  RightHadron=FS_RightHadron[sampledState];
706  } else
707  {
708  LeftHadron =FS_RightHadron[sampledState];
709  RightHadron=FS_LeftHadron[sampledState];
710  }
711  //G4cout<<"Selected "<<SampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
712  }
713  } else
714  {
715  if (string->DecayIsQuark() && string->StableIsQuark() )
716  { //... there are quarks on cluster ends
717 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
718  G4cout<<"Q Q string LastSplit"<<G4endl;
719 #endif
720  Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
721  } else
722  { //... there is a Diquark on one of the cluster ends
723 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
724  G4cout<<"DiQ Q string Last Split"<<G4endl;
725 #endif
726 
727  Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
728  }
729 
730  if(NumberOf_FS == 0) return false;
731  G4int sampledState = SampleState();
732  LeftHadron =FS_LeftHadron[sampledState];
733  RightHadron=FS_RightHadron[sampledState];
734 
735 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
736  G4cout<<"Selected LeftHad RightHad "<<sampledState<<" "
737  <<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
738 #endif
739 
740  } // End of if(!string->FourQuarkString())
741 
742  G4LorentzVector LeftMom, RightMom;
743  G4ThreeVector Pos;
744 
745  Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
746  &RightMom, RightHadron->GetPDGMass(),
747  StringMass);
748 
749  LeftMom.boost(ClusterVel);
750  RightMom.boost(ClusterVel);
751 
752  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
753  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
754 
755  return true;
756 
757 }
758 
759 //----------------------------------------------------------------------------------------------------------
760 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass,
761  G4LorentzVector* AntiMom, G4double AntiMass,
762  G4double InitialMass)
763 {
764  // ------ Sampling of momenta of 2 last produced hadrons --------------------
765  G4ThreeVector Pt;
766  G4double MassMt2, AntiMassMt2;
767  G4double AvailablePz, AvailablePz2;
768  G4double ProbIsotropy = sqr(sqr(938.0/InitialMass)); // Uzhi May 2015
769 #ifdef debug_LUNDfragmentation
770  G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
771  G4cout<<"Masses "<<InitialMass<<" "<<Mass<<" "<<AntiMass<<G4endl;
772 #endif
773 
774  if((Mass > 930. || AntiMass > 930.) && (G4UniformRand() < ProbIsotropy)) // Uzhi May 2015
775  { //If there is a baryon
776  // ----------------- Isotropic decay ------------------------------------
777  G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
778  sqr(2.*Mass*AntiMass);
779  G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
780  //G4cout<<"P for isotr decay "<<Pabs<<G4endl;
781 
782  //... sample unit vector
783  G4double pz =1. - 2.*G4UniformRand();
784  G4double st = std::sqrt(1. - pz * pz)*Pabs;
785  G4double phi = 2.*pi*G4UniformRand();
786  G4double px = st*std::cos(phi);
787  G4double py = st*std::sin(phi);
788  pz *= Pabs;
789 
790  Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
791  Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
792 
793  AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
794  AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
795  //G4int Uzhi; G4cin>>Uzhi;
796  }
797  else
798  {
799  const G4int maxNumberOfLoops = 1000;
800  G4int loopCounter = 0;
801  do
802  {
803  // GF 22-May-09, limit sampled pt to allowed range
804 
805  G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
806  G4double termab = 4*sqr(Mass*AntiMass);
807  G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
808  G4double pt2max=(termD*termD - termab )/ termN ;
809  //G4cout<<"Anis "<<pt2max<<" "<<(termD*termD-termab)/(4.*InitialMass*InitialMass)<<G4endl;
810 
811  Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
812  //G4cout<<"Sampl pt2 "<<Pt2<<G4endl;
813  MassMt2 = Mass * Mass + Pt2;
814  AntiMassMt2= AntiMass * AntiMass + Pt2;
815 
816  AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
817  4.*MassMt2*AntiMassMt2;
818  }
819  while( (AvailablePz2 < 0.) && // GF will occur only for numerical precision problem with limit in sampled pt
820  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
821 
822  if ( loopCounter >= maxNumberOfLoops ) {
823  AvailablePz2 = 0.0;
824  }
825 
826  AvailablePz2 /=(4.*InitialMass*InitialMass);
827  AvailablePz = std::sqrt(AvailablePz2);
828 
829  G4double Px=Pt.getX();
830  G4double Py=Pt.getY();
831 
832  Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
833  Mom->setE(std::sqrt(MassMt2+AvailablePz2));
834 
835  AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
836  AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
837  }
838 }
839 
840 //-----------------------------------------------------------------------------
841 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
842  G4FragmentingString * string,
843  G4FragmentingString * newString)
844 {
845  G4LorentzVector String4Momentum=string->Get4Momentum();
846  G4double StringMT2=string->MassT2();
847  G4double StringMT =std::sqrt(StringMT2);
848 
849  G4double HadronMass = pHadron->GetPDGMass();
850 
851  SetMinimalStringMass(newString);
852 
853 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
854  G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
855  G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()<<" "<<std::sqrt(StringMT2)<<G4endl;
856  G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
857  G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
858  <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
859 #endif
860 
861  if(HadronMass + MinimalStringMass > string->Mass())
862  {
863 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
864  G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
865 #endif
866  return 0;
867  }// have to start all over!
868 
869  String4Momentum.setPz(0.);
870  G4ThreeVector StringPt=String4Momentum.vect();
871 
872  // calculate and assign hadron transverse momentum component HadronPx and HadronPy
873  G4ThreeVector HadronPt , RemSysPt;
874  G4double HadronMassT2, ResidualMassT2;
875 
876  //... sample Pt of the hadron
877  G4int attempt=0;
878  do
879  {
880  attempt++; if(attempt > StringLoopInterrupt) return 0;
881 
882  HadronPt =SampleQuarkPt() + string->DecayPt();
883  HadronPt.setZ(0);
884  RemSysPt = StringPt - HadronPt;
885 
886  HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
887  ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
888 
889  } while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
890 
891  //... sample z to define hadron longitudinal momentum and energy
892  //... but first check the available phase space
893 
894  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
895  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
896 
897  if(Pz2 < 0 ) {return 0;} // have to start all over!
898 
899  //... then compute allowed z region z_min <= z <= z_max
900 
901  G4double Pz = std::sqrt(Pz2);
902  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
903 // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // Uzhi 2014
904  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
905 
906  if (zMin >= zMax) return 0; // have to start all over!
907 
908  G4double z = GetLightConeZ(zMin, zMax,
909  string->GetDecayParton()->GetPDGEncoding(), pHadron,
910  HadronPt.x(), HadronPt.y());
911 
912  //... now compute hadron longitudinal momentum and energy
913  // longitudinal hadron momentum component HadronPz
914 
915  HadronPt.setZ(0.5* string->GetDecayDirection() *
916  (z * string->LightConeDecay() -
917  HadronMassT2/(z * string->LightConeDecay())));
918  G4double HadronE = 0.5* (z * string->LightConeDecay() +
919  HadronMassT2/(z * string->LightConeDecay()));
920 
921  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
922 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
923  G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
924  G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
925  G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
926 //G4int Uzhi; G4cin>>Uzhi;
927  G4cout<<"Out of LUND SplitEandP "<<G4endl;
928 #endif
929 
930  return a4Momentum;
931 }
932 
933 //-----------------------------------------------------------------------------------------
934 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
935  G4int PDGEncodingOfDecayParton,
936  G4ParticleDefinition* pHadron,
937  G4double Px, G4double Py)
938 {
939  G4double Mass = pHadron->GetPDGMass();
940 // G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
941 
942  G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
943 
944  G4double alund;
945  G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
946  if(std::abs(PDGEncodingOfDecayParton) < 1000)
947  { // ---------------- Quark fragmentation ----------------------
948  alund=0.7/GeV/GeV;
949  // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
950  // const G4double blund = 1;
951 
952  zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
953  maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-alund*Mt2/zOfMaxyf);
954 
955  const G4int maxNumberOfLoops = 1000;
956  G4int loopCounter = 0;
957  do
958  {
959  z = zmin + G4UniformRand()*(zmax-zmin);
960  yf = (1-z)/z * G4Exp(-alund*Mt2/z);
961  // yf = G4Pow::GetInstance()->powA(1.0-z,blund)/z*G4Exp(-alund*Mt2/z
962  }
963  while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
964  if ( loopCounter >= maxNumberOfLoops ) {
965  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
966  }
967  return z;
968  }
969 
970  if(std::abs(PDGEncodingOfDecayParton) > 1000) // Uzhi Sept. 2014
971  {
972 /*
973  if(HadronEncoding < 3000)
974  {
975  maxYf=(zmax-zmin);
976  do
977  {
978  z = zmin + G4UniformRand()*(zmax-zmin);
979  //yf=G4Exp(-sqr(z-Zc)/2/sqr(0.28)); // 0.42 0.632 0.28 a'la UrQMD
980  yf =(z-zmin);
981  }
982  while (G4UniformRand()*maxYf > yf);
983  }
984  else
985  { // Strange baryons
986 */
987  z = zmin + G4UniformRand()*(zmax-zmin);
988 // }
989  }
990 
991  return z;
992 
993 }
994 
995 //------------------------------------------------------------------------
996 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
997 {
998  G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
999  return lam;
1000 }
1001 
1002 
1003 //------------------------------------------------------------------------
1004 //------------------------------------------------------------------------
1005 // Internal methods introduced to improve the code structure (AR Nov 2011)
1006 //------------------------------------------------------------------------
1007 //------------------------------------------------------------------------
1008 
1009 //----------------------------------------------------------------------------------------
1010 G4bool G4LundStringFragmentation::Loop_toFragmentString(G4ExcitedString * & theStringInCMS,
1011  G4KineticTrackVector * & LeftVector,
1012  G4KineticTrackVector * & RightVector)
1013 {
1014 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
1015  G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
1016  <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
1017  <<theStringInCMS->GetDirection()<< G4endl;
1018 #endif
1019 
1020  G4bool final_success=false;
1021  G4bool inner_success=true;
1022  G4int attempt=0;
1023  while ( ! final_success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
1024  { // If the string fragmentation do not be happend, repeat the fragmentation.
1025 
1026  G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
1027  //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
1028  // Cleaning up the previously produced hadrons
1029  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
1030  LeftVector->clear();
1031  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
1032  RightVector->clear();
1033 
1034  // Main fragmentation loop until the string will not be able to fragment
1035  inner_success=true; // set false on failure.
1036  const G4int maxNumberOfLoops = 1000;
1037  G4int loopCounter = -1;
1038  while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
1039  { // Split current string into hadron + new string
1040 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
1041  G4cout<<"The string can fragment. "<<G4endl;;
1042 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
1043 #endif
1044  G4FragmentingString *newString=0; // used as output from SplitUp.
1045  G4KineticTrack * Hadron=Splitup(currentString,newString);
1046  if ( Hadron != 0 ) // Store the hadron
1047  {
1048 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
1049  G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
1050 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
1051 #endif
1052 
1053  if ( currentString->GetDecayDirection() > 0 )
1054  {
1055  LeftVector->push_back(Hadron);
1056  } else
1057  {
1058  RightVector->push_back(Hadron);
1059  }
1060  delete currentString;
1061  currentString=newString;
1062  }
1063  };
1064  if ( loopCounter >= maxNumberOfLoops ) {
1065  inner_success=false;
1066  }
1067  // Split remaining string into 2 final hadrons.
1068 #ifdef debug_LUNDfragmentation // Uzhi 20.06.2014
1069  G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
1070 #endif
1071 // Uzhi Sept 2014
1072  if ( inner_success && SplitLast(currentString, LeftVector, RightVector) )
1073  {
1074  final_success=true;
1075  }
1076 //
1077 //final_success=true; // Uzhi Sept 2014
1078  delete currentString;
1079  } // End of the loop where we try to fragment the string.
1080  return final_success;
1081 }
1082 
1083 //----------------------------------------------------------------------------------------
1084 G4bool G4LundStringFragmentation::
1085 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
1086  G4ParticleDefinition * & LeftHadron,
1087  G4ParticleDefinition * & RightHadron)
1088 {
1089  G4double StringMass = string->Mass();
1090  G4int cClusterInterrupt = 0;
1091  do
1092  {
1093  //G4cout<<"cClusterInterrupt "<<cClusterInterrupt<<G4endl;
1094  if (cClusterInterrupt++ >= ClusterLoopInterrupt)
1095  {
1096  return false;
1097  }
1098 
1099  G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
1100  G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
1101 
1102  G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
1103  G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
1104 
1105  if(G4UniformRand()<0.5)
1106  {
1107  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
1108  FindParticle(RightQuark1));
1109  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
1110  FindParticle(RightQuark2));
1111  } else
1112  {
1113  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
1114  FindParticle(RightQuark2));
1115  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
1116  FindParticle(RightQuark1));
1117  }
1118 
1119  //... repeat procedure, if mass of cluster is too low to produce hadrons
1120  //... ClusterMassCut = 0.15*GeV model parameter
1121  }
1122  while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass())); /* Loop checking, 07.08.2015, A.Ribon */
1123 
1124  return true;
1125 }
1126 
1127 //----------------------------------------------------------------------------------------
1128 G4bool G4LundStringFragmentation::
1129 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
1130  G4ParticleDefinition * & LeftHadron,
1131  G4ParticleDefinition * & RightHadron)
1132 {
1133  // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
1134 
1135  G4double StringMass = string->Mass();
1136  G4double StringMassSqr= sqr(StringMass);
1137  G4ParticleDefinition * Di_Quark;
1138  G4ParticleDefinition * Anti_Di_Quark;
1139 
1140  if(string->GetLeftParton()->GetPDGEncoding() < 0)
1141  {
1142  Anti_Di_Quark =string->GetLeftParton();
1143  Di_Quark=string->GetRightParton();
1144  } else
1145  {
1146  Anti_Di_Quark =string->GetRightParton();
1147  Di_Quark=string->GetLeftParton();
1148  }
1149 
1150  G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
1151  G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
1152  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1153  G4int AbsIDdi_quark =std::abs(IDdi_quark);
1154 
1155  G4int ADi_q1=AbsIDAnti_di_quark/1000;
1156  G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
1157 
1158  G4int Di_q1=AbsIDdi_quark/1000;
1159  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1160 
1161  NumberOf_FS=0;
1162  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1163  {
1164  G4int StateADiQ=0;
1165  const G4int maxNumberOfLoops = 1000;
1166  G4int loopCounter = 0;
1167  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1168  {
1169  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
1170  -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
1171  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1172 
1173  //G4cout<<"Anti Bar "<<LeftHadron->GetParticleName()<<G4endl;
1174 
1175  G4int StateDiQ=0;
1176  const G4int maxNumberOfInternalLoops = 1000;
1177  G4int internalLoopCounter = 0;
1178  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1179  {
1180  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
1181  +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1182  G4double RightHadronMass=RightHadron->GetPDGMass();
1183 
1184  if(StringMass > LeftHadronMass + RightHadronMass)
1185  {
1186  if ( NumberOf_FS > 34 ) {
1187  G4ExceptionDescription ed;
1188  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1189  G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
1190  "HAD_LUND_001", JustWarning, ed );
1191  NumberOf_FS = 34;
1192  }
1193 
1194  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1195  sqr(RightHadronMass));
1196  //FS_Psqr=1.;
1197  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
1198  BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
1199  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1200  Prob_QQbar[ProdQ-1];
1201 
1202  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1203  FS_RightHadron[NumberOf_FS]= RightHadron;
1204 
1205  NumberOf_FS++;
1206  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
1207 
1208  StateDiQ++;
1209 
1210  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1211  ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1212  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1213  return false;
1214  }
1215 
1216  StateADiQ++;
1217  } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
1218  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1219  if ( loopCounter >= maxNumberOfLoops ) {
1220  return false;
1221  }
1222  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1223 
1224  return true;
1225 }
1226 
1227 //----------------------------------------------------------------------------------------
1228 G4bool G4LundStringFragmentation::
1229 Quark_Diquark_lastSplitting(G4FragmentingString * & string,
1230  G4ParticleDefinition * & LeftHadron,
1231  G4ParticleDefinition * & RightHadron)
1232 {
1233  G4double StringMass = string->Mass();
1234  G4double StringMassSqr= sqr(StringMass);
1235 
1236  G4ParticleDefinition * Di_Quark;
1237  G4ParticleDefinition * Quark;
1238 
1239  if(string->GetLeftParton()->GetParticleSubType()== "quark")
1240  {
1241  Quark =string->GetLeftParton();
1242  Di_Quark=string->GetRightParton();
1243  } else
1244  {
1245  Quark =string->GetRightParton();
1246  Di_Quark=string->GetLeftParton();
1247  }
1248 
1249  G4int IDquark =Quark->GetPDGEncoding();
1250  G4int AbsIDquark =std::abs(IDquark);
1251  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1252  G4int AbsIDdi_quark=std::abs(IDdi_quark);
1253  G4int Di_q1=AbsIDdi_quark/1000;
1254  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1255 
1256  G4int SignDiQ= 1;
1257  if(IDdi_quark < 0) SignDiQ=-1;
1258 
1259  NumberOf_FS=0;
1260  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1261  {
1262  G4int SignQ;
1263  if(IDquark > 0)
1264  { SignQ=-1;
1265  if(IDquark == 2) SignQ= 1;
1266  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1267  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1268  } else
1269  {
1270  SignQ= 1;
1271  if(IDquark == -2) SignQ=-1;
1272  if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1273  if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1274  }
1275 
1276  if(AbsIDquark == ProdQ) SignQ= 1;
1277 
1278  //G4cout<<G4endl;
1279  //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl;
1280 
1281  G4int StateQ=0;
1282  const G4int maxNumberOfLoops = 1000;
1283  G4int loopCounter = 0;
1284  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1285  {
1286  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1287  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1288  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1289 
1290  G4int StateDiQ=0;
1291  const G4int maxNumberOfInternalLoops = 1000;
1292  G4int internalLoopCounter = 0;
1293  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1294  {
1295  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1296  Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1297  G4double RightHadronMass=RightHadron->GetPDGMass();
1298 
1299  if(StringMass > LeftHadronMass + RightHadronMass)
1300  {
1301  if ( NumberOf_FS > 34 ) {
1302  G4ExceptionDescription ed;
1303  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1304  G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1305  "HAD_LUND_002", JustWarning, ed );
1306  NumberOf_FS = 34;
1307  }
1308 
1309  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1310  sqr(RightHadronMass));
1311  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1312  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1313  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1314  Prob_QQbar[ProdQ-1];
1315 
1316  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1317  FS_RightHadron[NumberOf_FS]= RightHadron;
1318 
1319  NumberOf_FS++;
1320  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
1321 
1322  StateDiQ++;
1323 
1324  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1325  ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1326  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1327  return false;
1328  }
1329 
1330  StateQ++;
1331  } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1332  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1333  if ( loopCounter >= maxNumberOfLoops ) {
1334  return false;
1335  }
1336  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1337 
1338  return true;
1339 }
1340 
1341 //----------------------------------------------------------------------------------------
1342 G4bool G4LundStringFragmentation::
1343 Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
1344  G4ParticleDefinition * & LeftHadron,
1345  G4ParticleDefinition * & RightHadron)
1346 {
1347  G4double StringMass = string->Mass();
1348  G4double StringMassSqr= sqr(StringMass);
1349 
1350  G4ParticleDefinition * Quark;
1351  G4ParticleDefinition * Anti_Quark;
1352 
1353  if(string->GetLeftParton()->GetPDGEncoding()>0)
1354  {
1355  Quark =string->GetLeftParton();
1356  Anti_Quark=string->GetRightParton();
1357  } else
1358  {
1359  Quark =string->GetRightParton();
1360  Anti_Quark=string->GetLeftParton();
1361  }
1362 
1363  G4int IDquark =Quark->GetPDGEncoding();
1364  G4int AbsIDquark =std::abs(IDquark);
1365  G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1366  G4int AbsIDanti_quark=std::abs(IDanti_quark);
1367 
1368  NumberOf_FS=0;
1369  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1370  {
1371  G4int SignQ=-1;
1372  if(IDquark == 2) SignQ= 1;
1373  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1374  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1375  if(IDquark == ProdQ) SignQ= 1;
1376 
1377  G4int SignAQ= 1;
1378  if(IDanti_quark == -2) SignAQ=-1;
1379  if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
1380  if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
1381  if(AbsIDanti_quark == ProdQ) SignAQ= 1;
1382 
1383  G4int StateQ=0;
1384  const G4int maxNumberOfLoops = 1000;
1385  G4int loopCounter = 0;
1386  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1387  {
1388  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1389  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1390  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1391 
1392  G4int StateAQ=0;
1393  const G4int maxNumberOfInternalLoops = 1000;
1394  G4int internalLoopCounter = 0;
1395  do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0);
1396  {
1397  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1398  Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1399  G4double RightHadronMass=RightHadron->GetPDGMass();
1400 
1401  if(StringMass > LeftHadronMass + RightHadronMass)
1402  {
1403  if ( NumberOf_FS > 34 ) {
1404  G4ExceptionDescription ed;
1405  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1406  G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1407  "HAD_LUND_003", JustWarning, ed );
1408  NumberOf_FS = 34;
1409  }
1410 
1411  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1412  sqr(RightHadronMass));
1413  //FS_Psqr=1.;
1414  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1415  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1416  MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1417  Prob_QQbar[ProdQ-1];
1418 
1419  if(string->GetLeftParton()->GetPDGEncoding()>0)
1420  {
1421  FS_LeftHadron[NumberOf_FS] = RightHadron;
1422  FS_RightHadron[NumberOf_FS]= LeftHadron;
1423  } else
1424  {
1425  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1426  FS_RightHadron[NumberOf_FS]= RightHadron;
1427  }
1428  NumberOf_FS++;
1429 
1430  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
1431 
1432  StateAQ++;
1433  } while( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1434  ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1435  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1436  return false;
1437  }
1438 
1439  StateQ++;
1440  } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1441  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1442  if ( loopCounter >= maxNumberOfLoops ) {
1443  return false;
1444  }
1445  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1446 
1447  return true;
1448 }
1449 
1450 //----------------------------------------------------------------------------------------------------------
1451 G4int G4LundStringFragmentation::SampleState(void)
1452 {
1453  if ( NumberOf_FS > 34 ) {
1454  G4ExceptionDescription ed;
1455  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1456  G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1457  NumberOf_FS = 34;
1458  }
1459 
1460  G4double SumWeights=0.;
1461 
1462  for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}// G4cout<<i<<" "<<FS_Weight[i]<<G4endl;}
1463 
1464  G4double ksi=G4UniformRand();
1465  G4double Sum=0.;
1466  G4int indexPosition = 0;
1467 
1468  for(G4int i=0; i<NumberOf_FS; i++)
1469  {
1470  Sum+=(FS_Weight[i]/SumWeights);
1471  indexPosition=i;
1472  if(Sum >= ksi) break;
1473  }
1474  return indexPosition;
1475 }
1476 // Uzhi June 2014 Insert from G4ExcitedStringDecay.cc
1477 //-----------------------------------------------------------------------------
1478 
1479 G4ParticleDefinition * G4LundStringFragmentation::DiQuarkSplitup(
1480  G4ParticleDefinition* decay,
1481  G4ParticleDefinition *&created)
1482 {
1483  //... can Diquark break or not?
1484  if (G4UniformRand() < DiquarkBreakProb ){
1485  //... Diquark break
1486 
1487  G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
1488  G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
1489  if (G4UniformRand() < 0.5)
1490  {
1491  G4int Swap = stableQuarkEncoding;
1492  stableQuarkEncoding = decayQuarkEncoding;
1493  decayQuarkEncoding = Swap;
1494  }
1495 
1496  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
1497  // if we have a quark, we need antiquark)
1498 
1499 //G4cout<<"GetStrangeSuppress() "<<GetStrangeSuppress()<<G4endl;
1500 //G4int Uzhi; G4cin>>Uzhi;
1501 
1502 //G4double StrSup=GetStrangeSuppress();
1503 //StrangeSuppress=0.34;
1504  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
1505 //StrangeSuppress=StrSup;
1506  //... Build new Diquark
1507  G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
1508  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1509  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1510  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
1511  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
1512  created = FindParticle(NewDecayEncoding);
1513  G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
1514  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
1515  return had;
1516 // return hadronizer->Build(QuarkPair.first, decayQuark);
1517 
1518  } else {
1519  //... Diquark does not break
1520 
1521  G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
1522  // if we have a diquark, we need quark)
1523 
1524 G4double StrSup=GetStrangeSuppress(); // Uzhi Sept. 2014
1525 StrangeSuppress=0.43; //0.42 0.38 // Uzhi Sept. 2014
1526  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
1527 StrangeSuppress=StrSup; // Uzhi Sept. 2014
1528 
1529  created = QuarkPair.second;
1530 
1531  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
1532  return had;
1533  }
1534 }
1535 // Uzhi June 2014 End of the inserting
G4LundStringFragmentation::FragmentString
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
Definition: G4LundStringFragmentation.cc:515
G4Exp.hh
G4FragmentingString::GetRightParton
G4ParticleDefinition * GetRightParton(void) const
Definition: G4FragmentingString.hh:143
G4FragmentingString.hh
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:566
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition: G4Parton.hh:124
C1
const double C1
Definition: MCGIDI_KalbachMann.cc:15
G4LundStringFragmentation::Quark_AntiQuark_lastSplitting
G4bool Quark_AntiQuark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:1343
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4LundStringFragmentation::SplitEandP
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
Definition: G4LundStringFragmentation.cc:841
G4LorentzRotation
CLHEP::HepLorentzRotation G4LorentzRotation
Definition: G4LorentzRotation.hh:35
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:230
G4VLongitudinalStringDecay::CPExcited
G4ExcitedString * CPExcited(const G4ExcitedString &string)
Definition: G4VLongitudinalStringDecay.cc:308
z
G4double z
Definition: TRTMaterials.hh:39
G4LundStringFragmentation::IsFragmentable
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:635
G4KineticTrack::SetFormationTime
void SetFormationTime(G4double aFormationTime)
Definition: G4KineticTrack.hh:225
G4VLongitudinalStringDecay::SetStrangenessSuppression
void SetStrangenessSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:556
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:140
G4HadronBuilder::Build
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition: G4HadronBuilder.cc:53
G4LundStringFragmentation::Quark_Diquark_lastSplitting
G4bool Quark_Diquark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:1229
G4ExcitedString::GetLeftParton
G4Parton * GetLeftParton(void) const
Definition: G4ExcitedString.cc:161
G4FragmentingString::Mass
G4double Mass() const
Definition: G4FragmentingString.cc:270
G4VLongitudinalStringDecay::SampleQuarkPt
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
Definition: G4VLongitudinalStringDecay.cc:499
G4VLongitudinalStringDecay::vectorMesonMix
std::vector< G4double > vectorMesonMix
Definition: G4VLongitudinalStringDecay.hh:189
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:181
G4VLongitudinalStringDecay::SetMassCut
virtual void SetMassCut(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:143
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4FragmentingString::GetDecayParton
G4ParticleDefinition * GetDecayParton() const
Definition: G4FragmentingString.hh:131
G4LundStringFragmentation::GetLightConeZ
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
Definition: G4LundStringFragmentation.cc:934
G4VLongitudinalStringDecay::CreatePartonPair
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
Definition: G4VLongitudinalStringDecay.cc:475
G4VLongitudinalStringDecay::SetStringTensionParameter
void SetStringTensionParameter(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:659
G4ExcitedString::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4ExcitedString.hh:149
G4LundStringFragmentation::Sample4Momentum
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
Definition: G4LundStringFragmentation.cc:760
G4FragmentingString::GetLeftParton
G4ParticleDefinition * GetLeftParton(void) const
Definition: G4FragmentingString.hh:137
G4VLongitudinalStringDecay::SetDiquarkBreakProbability
void SetDiquarkBreakProbability(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:578
G4LundStringFragmentation::SetMinimalStringMass2
void SetMinimalStringMass2(const G4double aValue)
Definition: G4LundStringFragmentation.cc:509
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4INCL::ClusterDecay::decay
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
Definition: G4INCLClusterDecay.cc:553
G4VLongitudinalStringDecay::LightFragmentationTest
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
Definition: G4VLongitudinalStringDecay.cc:149
G4ExcitedString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4ExcitedString.hh:161
G4LundStringFragmentation::Loop_toFragmentString
G4bool Loop_toFragmentString(G4ExcitedString *&theStringInCMS, G4KineticTrackVector *&LeftVector, G4KineticTrackVector *&RightVector)
Definition: G4LundStringFragmentation.cc:1010
G4KineticTrack::GetFormationTime
G4double GetFormationTime() const
Definition: G4KineticTrack.hh:220
G4VLongitudinalStringDecay::pDefPair
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
Definition: G4VLongitudinalStringDecay.hh:71
G4LundStringFragmentation::SetMinimalStringMass
void SetMinimalStringMass(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:415
G4bool
bool G4bool
Definition: G4Types.hh:79
G4ExcitedString::TransformToAlignedCms
G4LorentzRotation TransformToAlignedCms()
Definition: G4ExcitedString.hh:237
G4VLongitudinalStringDecay::FindParticle
G4ParticleDefinition * FindParticle(G4int Encoding)
Definition: G4VLongitudinalStringDecay.cc:281
G4KineticTrack::SetPosition
void SetPosition(const G4ThreeVector aPosition)
Definition: G4KineticTrack.hh:235
G4LundStringFragmentation::StopFragmenting
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:644
GeV
static const double GeV
Definition: G4SIunits.hh:214
Randomize.hh
G4FragmentingString::GetDecayDirection
G4int GetDecayDirection() const
Definition: G4FragmentingString.cc:190
G4KineticTrack::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4KineticTrack.hh:250
G4Quarks.hh
G4LundStringFragmentation.hh
G4LundStringFragmentation::lambda
G4double lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
Definition: G4LundStringFragmentation.cc:996
G4FragmentingString::FourQuarkString
G4bool FourQuarkString(void) const
Definition: G4FragmentingString.cc:200
G4VLongitudinalStringDecay::scalarMesonMix
std::vector< G4double > scalarMesonMix
Definition: G4VLongitudinalStringDecay.hh:190
G4LundStringFragmentation::DiQuarkSplitup
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
Definition: G4LundStringFragmentation.cc:1479
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4DiQuarks.hh
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ExcitedString::GetRightParton
G4Parton * GetRightParton(void) const
Definition: G4ExcitedString.cc:168
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting
G4bool Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:1129
pi
static const double pi
Definition: G4SIunits.hh:74
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4ExcitedString::GetTimeOfCreation
G4double GetTimeOfCreation() const
Definition: G4ExcitedString.hh:137
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4ExcitedString::GetDirection
G4int GetDirection(void) const
Definition: G4ExcitedString.cc:154
G4endl
#define G4endl
Definition: G4ios.hh:61
G4LundStringFragmentation::FS_LeftHadron
G4ParticleDefinition * FS_LeftHadron[35]
Definition: G4LundStringFragmentation.hh:133
G4LundStringFragmentation::Diquark_AntiDiquark_belowThreshold_lastSplitting
G4bool Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:1085
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4LundStringFragmentation::SplitLast
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
Definition: G4LundStringFragmentation.cc:659
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:240
G4Pow.hh
G4VLongitudinalStringDecay::CalculateHadronTimePosition
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
Definition: G4VLongitudinalStringDecay.cc:516
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:210
fermi
static const double fermi
Definition: G4SIunits.hh:102
G4VLongitudinalStringDecay::Splitup
G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
Definition: G4VLongitudinalStringDecay.cc:317
G4LundStringFragmentation::FS_RightHadron
G4ParticleDefinition * FS_RightHadron[35]
Definition: G4LundStringFragmentation.hh:133
G4VLongitudinalStringDecay::SetDiquarkSuppression
void SetDiquarkSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:567
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35