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