Geant4_10
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 76744 2013-11-14 15:50:20Z 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 // Class G4LundStringFragmentation
44 //*************************************************************************************
45 
46 G4LundStringFragmentation::G4LundStringFragmentation()
47 {
48 // ------ For estimation of a minimal string mass ---------------
49  Mass_of_light_quark =140.*MeV;
50  Mass_of_heavy_quark =500.*MeV;
51  Mass_of_string_junction=720.*MeV;
52 // ------ An estimated minimal string mass ----------------------
53  MinimalStringMass = 0.;
54  MinimalStringMass2 = 0.;
55 // ------ Minimal invariant mass used at a string fragmentation -
56  WminLUND = 0.45*GeV; //0.23*GeV; // Uzhi 0.7 -> 0.23 3.8.10 //0.8 1.5
57 // ------ Smooth parameter used at a string fragmentation for ---
58 // ------ smearinr sharp mass cut-off ---------------------------
59  SmoothParam = 0.2;
60 
61 // SetStringTensionParameter(0.25);
62  SetStringTensionParameter(1.);
63  SetDiquarkSuppression(0.087); // Uzhi 18.05.2012
64  SetDiquarkBreakProbability(0.05);
65  SetStrangenessSuppression(0.47); // Uzhi 18.05.2012
66 
67 // For treating of small string decays
68  for(G4int i=0; i<3; i++)
69  { for(G4int j=0; j<3; j++)
70  { for(G4int k=0; k<6; k++)
71  { Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
72  }
73  }
74  }
75 //--------------------------
76  Meson[0][0][0]=111; // dbar-d Pi0
77  MesonWeight[0][0][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
78 
79  Meson[0][0][1]=221; // dbar-d Eta
80  MesonWeight[0][0][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
81 
82  Meson[0][0][2]=331; // dbar-d EtaPrime
83  MesonWeight[0][0][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
84 
85  Meson[0][0][3]=113; // dbar-d Rho0
86  MesonWeight[0][0][3]=pspin_meson*(1.-vectorMesonMix[0]);
87 
88  Meson[0][0][4]=223; // dbar-d Omega
89  MesonWeight[0][0][4]=pspin_meson*(vectorMesonMix[0]);
90 //--------------------------
91 
92  Meson[0][1][0]=211; // dbar-u Pi+
93  MesonWeight[0][1][0]=(1.-pspin_meson);
94 
95  Meson[0][1][1]=213; // dbar-u Rho+
96  MesonWeight[0][1][1]=pspin_meson;
97 //--------------------------
98 
99  Meson[0][2][0]=311; // dbar-s K0bar
100  MesonWeight[0][2][0]=(1.-pspin_meson);
101 
102  Meson[0][2][1]=313; // dbar-s K*0bar
103  MesonWeight[0][2][1]=pspin_meson;
104 //--------------------------
105 //--------------------------
106  Meson[1][0][0]=211; // ubar-d Pi-
107  MesonWeight[1][0][0]=(1.-pspin_meson);
108 
109  Meson[1][0][1]=213; // ubar-d Rho-
110  MesonWeight[1][0][1]=pspin_meson;
111 //--------------------------
112 
113  Meson[1][1][0]=111; // ubar-u Pi0
114  MesonWeight[1][1][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
115 
116  Meson[1][1][1]=221; // ubar-u Eta
117  MesonWeight[1][1][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
118 
119  Meson[1][1][2]=331; // ubar-u EtaPrime
120  MesonWeight[1][1][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
121 
122  Meson[1][1][3]=113; // ubar-u Rho0
123  MesonWeight[1][1][3]=pspin_meson*(1.-vectorMesonMix[0]);
124 
125  Meson[1][1][4]=223; // ubar-u Omega
126  MesonWeight[1][1][4]=pspin_meson*(scalarMesonMix[0]);
127 //--------------------------
128 
129  Meson[1][2][0]=321; // ubar-s K-
130  MesonWeight[1][2][0]=(1.-pspin_meson);
131 
132  Meson[1][2][1]=323; // ubar-s K*-bar -
133  MesonWeight[1][2][1]=pspin_meson;
134 //--------------------------
135 //--------------------------
136 
137  Meson[2][0][0]=311; // sbar-d K0
138  MesonWeight[2][0][0]=(1.-pspin_meson);
139 
140  Meson[2][0][1]=313; // sbar-d K*0
141  MesonWeight[2][0][1]=pspin_meson;
142 //--------------------------
143 
144  Meson[2][1][0]=321; // sbar-u K+
145  MesonWeight[2][1][0]=(1.-pspin_meson);
146 
147  Meson[2][1][1]=323; // sbar-u K*+
148  MesonWeight[2][1][1]=pspin_meson;
149 //--------------------------
150 
151  Meson[2][2][0]=221; // sbar-s Eta
152  MesonWeight[2][2][0]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
153 
154  Meson[2][2][1]=331; // sbar-s EtaPrime
155  MesonWeight[2][2][1]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
156 
157  Meson[2][2][3]=333; // sbar-s EtaPrime
158  MesonWeight[2][2][3]=pspin_meson*(vectorMesonMix[5]);
159 //--------------------------
160 
161  for(G4int i=0; i<3; i++)
162  { for(G4int j=0; j<3; j++)
163  { for(G4int k=0; k<3; k++)
164  { for(G4int l=0; l<4; l++)
165  { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
166  }
167  }
168  }
169 
170  G4double pspin_barion_in=pspin_barion;
171  //pspin_barion=0.75;
172 //---------------------------------------
173  Baryon[0][0][0][0]=1114; // Delta-
174  BaryonWeight[0][0][0][0]=1.;
175 
176 //---------------------------------------
177  Baryon[0][0][1][0]=2112; // neutron
178  BaryonWeight[0][0][1][0]=1.-pspin_barion;
179 
180  Baryon[0][0][1][1]=2114; // Delta0
181  BaryonWeight[0][0][1][1]=pspin_barion;
182 
183 //---------------------------------------
184  Baryon[0][0][2][0]=3112; // Sigma-
185  BaryonWeight[0][0][2][0]=1.-pspin_barion;
186 
187  Baryon[0][0][2][1]=3114; // Sigma*-
188  BaryonWeight[0][0][2][1]=pspin_barion;
189 
190 //---------------------------------------
191  Baryon[0][1][0][0]=2112; // neutron
192  BaryonWeight[0][1][0][0]=1.-pspin_barion;
193 
194  Baryon[0][1][0][1]=2114; // Delta0
195  BaryonWeight[0][1][0][1]=pspin_barion;
196 
197 //---------------------------------------
198  Baryon[0][1][1][0]=2212; // proton
199  BaryonWeight[0][1][1][0]=1.-pspin_barion;
200 
201  Baryon[0][1][1][1]=2214; // Delta+
202  BaryonWeight[0][1][1][1]=pspin_barion;
203 
204 //---------------------------------------
205  Baryon[0][1][2][0]=3122; // Lambda
206  BaryonWeight[0][1][2][0]=(1.-pspin_barion)*0.5;
207 
208  Baryon[0][1][2][1]=3212; // Sigma0
209  BaryonWeight[0][1][2][1]=(1.-pspin_barion)*0.5;
210 
211  Baryon[0][1][2][2]=3214; // Sigma*0
212  BaryonWeight[0][1][2][2]=pspin_barion;
213 
214 //---------------------------------------
215  Baryon[0][2][0][0]=3112; // Sigma-
216  BaryonWeight[0][2][0][0]=1.-pspin_barion;
217 
218  Baryon[0][2][0][1]=3114; // Sigma*-
219  BaryonWeight[0][2][0][1]=pspin_barion;
220 
221 //---------------------------------------
222  Baryon[0][2][1][0]=3122; // Lambda
223  BaryonWeight[0][2][1][0]=(1.-pspin_barion)*0.5;
224 
225  Baryon[0][2][1][1]=3212; // Sigma0
226  BaryonWeight[0][2][1][1]=(1.-pspin_barion)*0.5;
227 
228  Baryon[0][2][1][2]=3214; // Sigma*0
229  BaryonWeight[0][2][1][2]=pspin_barion;
230 
231 //---------------------------------------
232  Baryon[0][2][2][0]=3312; // Theta-
233  BaryonWeight[0][2][2][0]=1.-pspin_barion;
234 
235  Baryon[0][2][2][1]=3314; // Theta*-
236  BaryonWeight[0][2][2][1]=pspin_barion;
237 
238 //---------------------------------------
239 //---------------------------------------
240  Baryon[1][0][0][0]=2112; // neutron
241  BaryonWeight[1][0][0][0]=1.-pspin_barion;
242 
243  Baryon[1][0][0][1]=2114; // Delta0
244  BaryonWeight[1][0][0][1]=pspin_barion;
245 
246 //---------------------------------------
247  Baryon[1][0][1][0]=2212; // proton
248  BaryonWeight[1][0][1][0]=1.-pspin_barion;
249 
250  Baryon[1][0][1][1]=2214; // Delta+
251  BaryonWeight[1][0][1][1]=pspin_barion;
252 
253 //---------------------------------------
254  Baryon[1][0][2][0]=3122; // Lambda
255  BaryonWeight[1][0][2][0]=(1.-pspin_barion)*0.5;
256 
257  Baryon[1][0][2][1]=3212; // Sigma0
258  BaryonWeight[1][0][2][1]=(1.-pspin_barion)*0.5;
259 
260  Baryon[1][0][2][2]=3214; // Sigma*0
261  BaryonWeight[1][0][2][2]=pspin_barion;
262 
263 //---------------------------------------
264  Baryon[1][1][0][0]=2212; // proton
265  BaryonWeight[1][1][0][0]=1.-pspin_barion;
266 
267  Baryon[1][1][0][1]=2214; // Delta+
268  BaryonWeight[1][1][0][1]=pspin_barion;
269 
270 //---------------------------------------
271  Baryon[1][1][1][0]=2224; // Delta++
272  BaryonWeight[1][1][1][0]=1.;
273 
274 //---------------------------------------
275  Baryon[1][1][2][0]=3222; // Sigma+
276  BaryonWeight[1][1][2][0]=1.-pspin_barion;
277 
278  Baryon[1][1][2][1]=3224; // Sigma*+
279  BaryonWeight[1][1][2][1]=pspin_barion;
280 
281 //---------------------------------------
282  Baryon[1][2][0][0]=3122; // Lambda
283  BaryonWeight[1][2][0][0]=(1.-pspin_barion)*0.5;
284 
285  Baryon[1][2][0][1]=3212; // Sigma0
286  BaryonWeight[1][2][0][1]=(1.-pspin_barion)*0.5;
287 
288  Baryon[1][2][0][2]=3214; // Sigma*0
289  BaryonWeight[1][2][0][2]=pspin_barion;
290 
291 //---------------------------------------
292  Baryon[1][2][1][0]=3222; // Sigma+
293  BaryonWeight[1][2][1][0]=1.-pspin_barion;
294 
295  Baryon[1][2][1][1]=3224; // Sigma*+
296  BaryonWeight[1][2][1][1]=pspin_barion;
297 
298 //---------------------------------------
299  Baryon[1][2][2][0]=3322; // Theta0
300  BaryonWeight[1][2][2][0]=1.-pspin_barion;
301 
302  Baryon[1][2][2][1]=3324; // Theta*0
303  BaryonWeight[1][2][2][1]=pspin_barion;
304 
305 //---------------------------------------
306 //---------------------------------------
307  Baryon[2][0][0][0]=3112; // Sigma-
308  BaryonWeight[2][0][0][0]=1.-pspin_barion;
309 
310  Baryon[2][0][0][1]=3114; // Sigma*-
311  BaryonWeight[2][0][0][1]=pspin_barion;
312 
313 //---------------------------------------
314  Baryon[2][0][1][0]=3122; // Lambda
315  BaryonWeight[2][0][1][0]=(1.-pspin_barion)*0.5;
316 
317  Baryon[2][0][1][1]=3212; // Sigma0
318  BaryonWeight[2][0][1][1]=(1.-pspin_barion)*0.5;
319 
320  Baryon[2][0][1][2]=3214; // Sigma*0
321  BaryonWeight[2][0][1][2]=pspin_barion;
322 
323 //---------------------------------------
324  Baryon[2][0][2][0]=3312; // Sigma-
325  BaryonWeight[2][0][2][0]=1.-pspin_barion;
326 
327  Baryon[2][0][2][1]=3314; // Sigma*-
328  BaryonWeight[2][0][2][1]=pspin_barion;
329 
330 //---------------------------------------
331  Baryon[2][1][0][0]=3122; // Lambda
332  BaryonWeight[2][1][0][0]=(1.-pspin_barion)*0.5;
333 
334  Baryon[2][1][0][1]=3212; // Sigma0
335  BaryonWeight[2][1][0][1]=(1.-pspin_barion)*0.5;
336 
337  Baryon[2][1][0][2]=3214; // Sigma*0
338  BaryonWeight[2][1][0][2]=pspin_barion;
339 
340 //---------------------------------------
341  Baryon[2][1][1][0]=3222; // Sigma+
342  BaryonWeight[2][1][1][0]=1.-pspin_barion;
343 
344  Baryon[2][1][1][1]=3224; // Sigma*+
345  BaryonWeight[2][1][1][1]=pspin_barion;
346 
347 //---------------------------------------
348  Baryon[2][1][2][0]=3322; // Theta0
349  BaryonWeight[2][1][2][0]=1.-pspin_barion;
350 
351  Baryon[2][1][2][1]=3324; // Theta*0
352  BaryonWeight[2][1][2][1]=pspin_barion;
353 
354 //---------------------------------------
355  Baryon[2][2][0][0]=3312; // Theta-
356  BaryonWeight[2][2][0][0]=1.-pspin_barion;
357 
358  Baryon[2][2][0][1]=3314; // Theta*-
359  BaryonWeight[2][2][0][1]=pspin_barion;
360 
361 //---------------------------------------
362  Baryon[2][2][1][0]=3322; // Theta0
363  BaryonWeight[2][2][1][0]=1.-pspin_barion;
364 
365  Baryon[2][2][1][1]=3324; // Theta*0
366  BaryonWeight[2][2][1][1]=pspin_barion;
367 
368 //---------------------------------------
369  Baryon[2][2][2][0]=3334; // Omega
370  BaryonWeight[2][2][2][0]=1.;
371 
372 //---------------------------------------
373  pspin_barion=pspin_barion_in;
374  /*
375  for(G4int i=0; i<3; i++)
376  { for(G4int j=0; j<3; j++)
377  { for(G4int k=0; k<3; k++)
378  { for(G4int l=0; l<4; l++)
379  { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;}
380  }
381  }
382  }
383  G4int Uzhi;
384  G4cin>>Uzhi;
385  */
386  //StrangeSuppress=0.38;
387  Prob_QQbar[0]=StrangeSuppress; // Probability of ddbar production
388  Prob_QQbar[1]=StrangeSuppress; // Probability of uubar production
389  Prob_QQbar[2]=StrangeSuppress/(2.+StrangeSuppress);//(1.-2.*StrangeSuppress); // Probability of ssbar production
390 
391  //A.R. 25-Jul-2012 : Coverity fix.
392  for ( G4int i=0 ; i<35 ; i++ ) {
393  FS_LeftHadron[i] = 0;
394  FS_RightHadron[i] = 0;
395  FS_Weight[i] = 0.0;
396  }
397  NumberOf_FS = 0;
398 
399 }
400 
401 // --------------------------------------------------------------
402 G4LundStringFragmentation::~G4LundStringFragmentation()
403 {}
404 
405 
406 //--------------------------------------------------------------------------------------
407 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
408 {
409  G4double EstimatedMass=0.;
410  G4int Number_of_quarks=0;
411 
412  G4double StringM=string->Get4Momentum().mag();
413 
414  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
415 
416  //G4cout<<"Min mass Qleft -------------------"<<G4endl;
417  //G4cout<<"String mass"<<string->Get4Momentum().mag()<<G4endl;
418  if( Qleft > 1000)
419  {
420  Number_of_quarks+=2;
421  G4int q1=Qleft/1000;
422  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
423  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
424 
425  G4int q2=(Qleft/100)%10;
426  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
427  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
428  EstimatedMass +=Mass_of_string_junction;
429  }
430  else
431  {
432  Number_of_quarks++;
433  if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
434  if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;}
435  }
436 
437  //G4cout<<"Min mass Qleft "<<Qleft<<" "<<EstimatedMass<<G4endl;
438 
439  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
440 
441  if( Qright > 1000)
442  {
443  Number_of_quarks+=2;
444  G4int q1=Qright/1000;
445  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
446  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
447 
448  G4int q2=(Qright/100)%10;
449  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
450  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
451  EstimatedMass +=Mass_of_string_junction;
452  }
453  else
454  {
455  Number_of_quarks++;
456  if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
457  if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;}
458  }
459 
460  //G4cout<<"Min mass Qright "<<Qright<<" "<<EstimatedMass<<G4endl;
461 
462  if(Number_of_quarks==2){EstimatedMass +=100.*MeV;}
463  if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
464  if(Number_of_quarks==4)
465  {
466  if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2020.;}//1880.;}
467  // if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2051.;}
468  else if((StringM > 2232.) && ( EstimatedMass < 2730)){EstimatedMass = 2570.;}
469  else if((StringM > 5130.) && ( EstimatedMass < 3450)){EstimatedMass = 5130.;}
470  else
471  {
472  EstimatedMass -=2.*Mass_of_string_junction;
473  if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
474  else {EstimatedMass+=100.*MeV;}
475  }
476  }
477 
478  //G4cout<<"EstimatedMass "<<EstimatedMass <<G4endl;
479  //G4int Uzhi; G4cin>>Uzhi;
480  MinimalStringMass=EstimatedMass;
481  SetMinimalStringMass2(EstimatedMass);
482 }
483 
484 //--------------------------------------------------------------------------------------
485 void G4LundStringFragmentation::SetMinimalStringMass2(
486  const G4double aValue)
487 {
488  MinimalStringMass2=aValue * aValue;
489 }
490 
491 //--------------------------------------------------------------------------------------
492 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(
493  const G4ExcitedString& theString)
494 {
495  // Can no longer modify Parameters for Fragmentation.
496  PastInitPhase=true;
497 
498  SetMassCut(160.*MeV); // For LightFragmentationTest it is required
499  // that no one pi-meson can be produced.
500 
501  G4FragmentingString aString(theString);
502  SetMinimalStringMass(&aString);
503 
504  G4KineticTrackVector * LeftVector(0);
505 
506  if(!IsFragmentable(&aString)) // produce 1 hadron
507  {
508  //G4cout<<"Non fragmentable"<<G4endl;
509  SetMassCut(1000.*MeV);
510  LeftVector=LightFragmentationTest(&theString);
511  SetMassCut(160.*MeV);
512  } // end of if(!IsFragmentable(&aString))
513 
514  if ( LeftVector != 0 ) {
515  // Uzhi insert 6.05.08 start
516  LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
517  LeftVector->operator[](0)->SetPosition(theString.GetPosition());
518  if(LeftVector->size() > 1)
519  {
520  // 2 hadrons created from qq-qqbar are stored
521  LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
522  LeftVector->operator[](1)->SetPosition(theString.GetPosition());
523  }
524  return LeftVector;
525  }
526 
527  // The string can fragment. At least two particles can be produced.
528  LeftVector =new G4KineticTrackVector;
529  G4KineticTrackVector * RightVector=new G4KineticTrackVector;
530 
531  G4ExcitedString *theStringInCMS=CPExcited(theString);
532  G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
533 
534  G4bool success = Loop_toFragmentString(theStringInCMS, LeftVector, RightVector);
535 
536  delete theStringInCMS;
537 
538  if ( ! success )
539  {
540  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
541  LeftVector->clear();
542  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
543  delete RightVector;
544  return LeftVector;
545  }
546 
547  // Join Left- and RightVector into LeftVector in correct order.
548  while(!RightVector->empty())
549  {
550  LeftVector->push_back(RightVector->back());
551  RightVector->erase(RightVector->end()-1);
552  }
553  delete RightVector;
554 
555  CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
556 
557  G4LorentzRotation toObserverFrame(toCms.inverse());
558 
559  G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
560  G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
561 
562  //G4cout<<"# prod hadrons "<<LeftVector->size()<<G4endl;
563  for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
564  {
565  G4KineticTrack* Hadron = LeftVector->operator[](C1);
566  G4LorentzVector Momentum = Hadron->Get4Momentum();
567  //G4cout<<"Hadron "<<Hadron->GetDefinition()->GetParticleName()<<" "<<Momentum<<G4endl;
568  Momentum = toObserverFrame*Momentum;
569  Hadron->Set4Momentum(Momentum);
570 
571  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
572  Momentum = toObserverFrame*Coordinate;
573  Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
574  G4ThreeVector aPosition(Momentum.vect());
575  Hadron->SetPosition(PositionOftheStringCreation+aPosition);
576  };
577 
578  return LeftVector;
579 }
580 
581 //----------------------------------------------------------------------------------
582 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
583 {
584  SetMinimalStringMass(string);
585  // return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
586  return MinimalStringMass < string->Get4Momentum().mag(); // 21.07.2010
587 }
588 
589 //----------------------------------------------------------------------------------------
590 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
591 {
592  SetMinimalStringMass(string);
593 
594  if (string->FourQuarkString())
595  {
596  return G4UniformRand() < std::exp(-0.0005*(string->Mass() - MinimalStringMass));
597  } else {
598  return G4UniformRand() < std::exp(-0.88e-6*(string->Mass()*string->Mass() -
599  MinimalStringMass*MinimalStringMass));
600  }
601 }
602 
603 //----------------------------------------------------------------------------------------------------------
604 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
605  G4KineticTrackVector * LeftVector,
606  G4KineticTrackVector * RightVector)
607 {
608  //... perform last cluster decay
609  //G4cout<<"Split last-----------------------------------------"<<G4endl;
610  G4LorentzVector Str4Mom=string->Get4Momentum();
611  G4ThreeVector ClusterVel=string->Get4Momentum().boostVector();
612  G4double StringMass=string->Mass();
613 
614  G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
615 
616  NumberOf_FS=0;
617  for(G4int i=0; i<35; i++) {FS_Weight[i]=0.;}
618 
619  //G4cout<<"StrMass "<<StringMass<<" q "<<string->GetLeftParton()->GetParticleName()<<" "<<string->GetRightParton()->GetParticleName()<<" StringMassSqr "<<StringMassSqr<<G4endl;
620 
621  string->SetLeftPartonStable(); // to query quark contents..
622 
623  if (string->FourQuarkString() )
624  {
625  // The string is qq-qqbar type. Diquarks are on the string ends
626  //G4cout<<"The string is qq-qqbar type. Diquarks are on the string ends"<<G4endl;
627 
628  if(StringMass-MinimalStringMass < 0.)
629  {
630  if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
631  {
632  return false;
633  }
634  } else
635  {
636  Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
637 
638  if(NumberOf_FS == 0) return false;
639  //G4cout<<"NumberOf_FS "<<NumberOf_FS<<G4endl;
640  G4int sampledState = SampleState();
641  //SampledState=16;
642  if(string->GetLeftParton()->GetPDGEncoding() < 0)
643  {
644  LeftHadron =FS_LeftHadron[sampledState];
645  RightHadron=FS_RightHadron[sampledState];
646  } else
647  {
648  LeftHadron =FS_RightHadron[sampledState];
649  RightHadron=FS_LeftHadron[sampledState];
650  }
651  //G4cout<<"Selected "<<SampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
652  }
653  } else
654  {
655  if (string->DecayIsQuark() && string->StableIsQuark() )
656  { //... there are quarks on cluster ends
657  //G4cout<<"Q Q string"<<G4endl;
658  Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
659  } else
660  { //... there is a Diquark on one of the cluster ends
661  //G4cout<<"DiQ Q string"<<G4endl;
662  Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
663  }
664 
665  if(NumberOf_FS == 0) return false;
666  //G4cout<<"NumberOf_FS "<<NumberOf_FS<<G4endl;
667  G4int sampledState = SampleState();
668  //sampledState=17;
669  LeftHadron =FS_LeftHadron[sampledState];
670  RightHadron=FS_RightHadron[sampledState];
671  //G4cout<<"Selected "<<sampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
672  //G4int Uzhi; G4cin>>Uzhi;
673 
674  } // End of if(!string->FourQuarkString())
675 
676  G4LorentzVector LeftMom, RightMom;
677  G4ThreeVector Pos;
678 
679  Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
680  &RightMom, RightHadron->GetPDGMass(),
681  StringMass);
682 
683  LeftMom.boost(ClusterVel);
684  RightMom.boost(ClusterVel);
685 
686  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
687  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
688 
689  return true;
690 
691 }
692 
693 //----------------------------------------------------------------------------------------------------------
694 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
695 {
696  // ------ Sampling of momenta of 2 last produced hadrons --------------------
697  G4ThreeVector Pt;
698  G4double MassMt2, AntiMassMt2;
699  G4double AvailablePz, AvailablePz2;
700 
701  //G4cout<<"Masses "<<InitialMass<<" "<<Mass<<" "<<AntiMass<<G4endl;
702  //
703 
704  if((Mass > 930. || AntiMass > 930.)) //If there is a baryon
705  {
706  // ----------------- Isotropic decay ------------------------------------
707  G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
708  sqr(2.*Mass*AntiMass);
709  G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
710  //G4cout<<"P for isotr decay "<<Pabs<<G4endl;
711 
712  //... sample unit vector
713  G4double pz =1. - 2.*G4UniformRand();
714  G4double st = std::sqrt(1. - pz * pz)*Pabs;
715  G4double phi = 2.*pi*G4UniformRand();
716  G4double px = st*std::cos(phi);
717  G4double py = st*std::sin(phi);
718  pz *= Pabs;
719 
720  Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
721  Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
722 
723  AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
724  AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
725  //G4int Uzhi; G4cin>>Uzhi;
726  }
727  else
728  //
729  {
730  do
731  {
732  // GF 22-May-09, limit sampled pt to allowed range
733 
734  G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
735  G4double termab = 4*sqr(Mass*AntiMass);
736  G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
737  G4double pt2max=(termD*termD - termab )/ termN ;
738  //G4cout<<"Anis "<<pt2max<<" "<<(termD*termD-termab)/(4.*InitialMass*InitialMass)<<G4endl;
739 
740  Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
741  //G4cout<<"Sampl pt2 "<<Pt2<<G4endl;
742  MassMt2 = Mass * Mass + Pt2;
743  AntiMassMt2= AntiMass * AntiMass + Pt2;
744 
745  AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
746  4.*MassMt2*AntiMassMt2;
747  }
748  while(AvailablePz2 < 0.); // GF will occur only for numerical precision problem with limit in sampled pt
749 
750  AvailablePz2 /=(4.*InitialMass*InitialMass);
751  AvailablePz = std::sqrt(AvailablePz2);
752 
753  G4double Px=Pt.getX();
754  G4double Py=Pt.getY();
755 
756  Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
757  Mom->setE(std::sqrt(MassMt2+AvailablePz2));
758 
759  AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
760  AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
761  }
762 }
763 
764 //-----------------------------------------------------------------------------
765 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
766  G4FragmentingString * string, G4FragmentingString * newString)
767 {
768  //G4cout<<"Start SplitEandP "<<G4endl;
769  G4LorentzVector String4Momentum=string->Get4Momentum();
770  G4double StringMT2=string->Get4Momentum().mt2();
771 
772  G4double HadronMass = pHadron->GetPDGMass();
773 
774  SetMinimalStringMass(newString);
775  //G4cout<<"HadM MinimalStringMassLeft StringM "<<HadronMass<<" "<<MinimalStringMass<<" "<<String4Momentum.mag()<<G4endl;
776 
777  if(HadronMass + MinimalStringMass > String4Momentum.mag()) {return 0;}// have to start all over!
778  String4Momentum.setPz(0.);
779  G4ThreeVector StringPt=String4Momentum.vect();
780 
781  // calculate and assign hadron transverse momentum component HadronPx and HadronPy
782  G4ThreeVector thePt;
783  thePt=SampleQuarkPt();
784 
785  G4ThreeVector HadronPt = thePt +string->DecayPt();
786  HadronPt.setZ(0);
787 
788  G4ThreeVector RemSysPt = StringPt - HadronPt;
789 
790  //... sample z to define hadron longitudinal momentum and energy
791  //... but first check the available phase space
792 
793  G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
794  G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
795 
796  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
797  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
798  //G4cout<<"Pz2 "<<Pz2<<G4endl;
799  if(Pz2 < 0 ) {return 0;} // have to start all over!
800 
801  //... then compute allowed z region z_min <= z <= z_max
802 
803  G4double Pz = std::sqrt(Pz2);
804  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
805  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
806 
807  //G4cout<<"if (zMin >= zMax) return 0 "<<zMin<<" "<<zMax<<G4endl;
808  if (zMin >= zMax) return 0; // have to start all over!
809 
810  G4double z = GetLightConeZ(zMin, zMax,
811  string->GetDecayParton()->GetPDGEncoding(), pHadron,
812  HadronPt.x(), HadronPt.y());
813  //G4cout<<"z "<<z<<G4endl;
814  //... now compute hadron longitudinal momentum and energy
815  // longitudinal hadron momentum component HadronPz
816 
817  HadronPt.setZ(0.5* string->GetDecayDirection() *
818  (z * string->LightConeDecay() -
819  HadronMassT2/(z * string->LightConeDecay())));
820 
821  G4double HadronE = 0.5* (z * string->LightConeDecay() +
822  HadronMassT2/(z * string->LightConeDecay()));
823 
824  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
825  //G4cout<<"Out SplitEandP "<<G4endl;
826  return a4Momentum;
827 }
828 
829 //-----------------------------------------------------------------------------------------
830 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
831  G4int PDGEncodingOfDecayParton,
832  G4ParticleDefinition* pHadron,
833  G4double Px, G4double Py)
834 {
835 
836  // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
837  // const G4double blund = 1;
838 
839  G4double Mass = pHadron->GetPDGMass();
840  // G4int HadronEncoding=pHadron->GetPDGEncoding();
841 
842  G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
843 
844  G4double alund;
845  if(std::abs(PDGEncodingOfDecayParton) < 1000)
846  { // ---------------- Quark fragmentation ----------------------
847  alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered
848  }
849  else
850  { // ---------------- Di-quark fragmentation ----------------------
851  alund=0.7/GeV/GeV; // 0.7 2.0
852  }
853  G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
854  G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
855  G4double z, yf;
856  do
857  {
858  z = zmin + G4UniformRand()*(zmax-zmin);
859  // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
860  yf = (1-z)/z * std::exp(-alund*Mt2/z);
861  }
862  while (G4UniformRand()*maxYf > yf);
863 
864 
865  return z;
866 }
867 
868 //------------------------------------------------------------------------
869 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
870 {
871  G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
872  return lam;
873 }
874 
875 
876 //------------------------------------------------------------------------
877 //------------------------------------------------------------------------
878 // Internal methods introduced to improve the code structure (AR Nov 2011)
879 //------------------------------------------------------------------------
880 //------------------------------------------------------------------------
881 
882 //----------------------------------------------------------------------------------------
883 G4bool G4LundStringFragmentation::Loop_toFragmentString(G4ExcitedString * & theStringInCMS,
884  G4KineticTrackVector * & LeftVector,
885  G4KineticTrackVector * & RightVector)
886 {
887  G4bool final_success=false;
888  G4bool inner_success=true;
889  G4int attempt=0;
890  while ( ! final_success && attempt++ < StringLoopInterrupt )
891  { // If the string fragmentation do not be happend, repeat the fragmentation.
892  G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
893  //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
894  // Cleaning up the previously produced hadrons
895  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
896  LeftVector->clear();
897  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
898  RightVector->clear();
899 
900  // Main fragmentation loop until the string will not be able to fragment
901  inner_success=true; // set false on failure.
902  while (! StopFragmenting(currentString) )
903  { // Split current string into hadron + new string
904  G4FragmentingString *newString=0; // used as output from SplitUp.
905  G4KineticTrack * Hadron=Splitup(currentString,newString);
906  if ( Hadron != 0 ) // Store the hadron
907  {
908  //G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
909  if ( currentString->GetDecayDirection() > 0 )
910  {
911  LeftVector->push_back(Hadron);
912  } else
913  {
914  RightVector->push_back(Hadron);
915  }
916  delete currentString;
917  currentString=newString;
918  }
919  };
920  // Split remaining string into 2 final hadrons.
921  if ( inner_success && SplitLast(currentString, LeftVector, RightVector) )
922  {
923  final_success=true;
924  }
925  delete currentString;
926  } // End of the loop where we try to fragment the string.
927  return final_success;
928 }
929 
930 //----------------------------------------------------------------------------------------
931 G4bool G4LundStringFragmentation::
932 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
933  G4ParticleDefinition * & LeftHadron,
934  G4ParticleDefinition * & RightHadron)
935 {
936  G4double StringMass = string->Mass();
937  G4int cClusterInterrupt = 0;
938  do
939  {
940  //G4cout<<"cClusterInterrupt "<<cClusterInterrupt<<G4endl;
941  if (cClusterInterrupt++ >= ClusterLoopInterrupt)
942  {
943  return false;
944  }
945 
946  G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
947  G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
948 
949  G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
950  G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
951 
952  if(G4UniformRand()<0.5)
953  {
954  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
955  FindParticle(RightQuark1));
956  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
957  FindParticle(RightQuark2));
958  } else
959  {
960  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
961  FindParticle(RightQuark2));
962  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
963  FindParticle(RightQuark1));
964  }
965 
966  //... repeat procedure, if mass of cluster is too low to produce hadrons
967  //... ClusterMassCut = 0.15*GeV model parameter
968  }
969  while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()));
970 
971  return true;
972 }
973 
974 //----------------------------------------------------------------------------------------
975 G4bool G4LundStringFragmentation::
976 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
977  G4ParticleDefinition * & LeftHadron,
978  G4ParticleDefinition * & RightHadron)
979 {
980  // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
981 
982  //G4cout<<"DiQ Anti-DiQ string"<<G4endl;
983  //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetRightParton()->GetPDGEncoding()<<G4endl;
984 
985  G4double StringMass = string->Mass();
986  G4double StringMassSqr= sqr(StringMass);
987  G4ParticleDefinition * Di_Quark;
988  G4ParticleDefinition * Anti_Di_Quark;
989 
990  if(string->GetLeftParton()->GetPDGEncoding() < 0)
991  {
992  Anti_Di_Quark =string->GetLeftParton();
993  Di_Quark=string->GetRightParton();
994  } else
995  {
996  Anti_Di_Quark =string->GetRightParton();
997  Di_Quark=string->GetLeftParton();
998  }
999 
1000  G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
1001  G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
1002  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1003  G4int AbsIDdi_quark =std::abs(IDdi_quark);
1004 
1005  G4int ADi_q1=AbsIDAnti_di_quark/1000;
1006  G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
1007 
1008  G4int Di_q1=AbsIDdi_quark/1000;
1009  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1010 
1011  //G4cout<<"IDs Anti "<<ADi_q1<<" "<<ADi_q2<<G4endl;
1012  //G4cout<<"IDs "<< Di_q1<<" "<< Di_q2<<G4endl;
1013 
1014  NumberOf_FS=0;
1015  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1016  {
1017  //G4cout<<"Insert QQbar "<<ProdQ<<G4endl;
1018 
1019  G4int StateADiQ=0;
1020  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1021  {
1022  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1023  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1024  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
1025  -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
1026  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1027 
1028  //G4cout<<"Anti Bar "<<LeftHadron->GetParticleName()<<G4endl;
1029 
1030  G4int StateDiQ=0;
1031  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1032  {
1033  //G4cout<<G4endl<<"Di_q1 Di_q2 ProdQ StateDiQ "<<BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1034  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
1035  +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1036  G4double RightHadronMass=RightHadron->GetPDGMass();
1037 
1038  //G4cout<<"Baryon "<<RightHadron->GetParticleName()<<G4endl;
1039 
1040  //G4cout<<"StringMass LeftHadronMass RightHadronMass "<<StringMass<<" "<<LeftHadronMass<<" "<< RightHadronMass<<G4endl;
1041 
1042  if(StringMass >= LeftHadronMass + RightHadronMass)
1043  {
1044  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1045  sqr(RightHadronMass));
1046  //FS_Psqr=1.;
1047  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
1048  BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
1049  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1050  Prob_QQbar[ProdQ-1];
1051 
1052  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1053  FS_RightHadron[NumberOf_FS]= RightHadron;
1054 
1055  //G4cout<<"State "<<NumberOf_FS<<" "<<BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<G4endl;
1056  //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl;
1057  NumberOf_FS++;
1058 
1059  if(NumberOf_FS > 34)
1060  {G4int Uzhi; G4cout<<"QQ_QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;}
1061  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1062 
1063  StateDiQ++;
1064  //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1065  } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0);
1066 
1067  StateADiQ++;
1068  } while(Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0);
1069  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1070 
1071  return true;
1072 }
1073 
1074 //----------------------------------------------------------------------------------------
1075 G4bool G4LundStringFragmentation::
1076 Quark_Diquark_lastSplitting(G4FragmentingString * & string,
1077  G4ParticleDefinition * & LeftHadron,
1078  G4ParticleDefinition * & RightHadron)
1079 {
1080  G4double StringMass = string->Mass();
1081  G4double StringMassSqr= sqr(StringMass);
1082 
1083  G4ParticleDefinition * Di_Quark;
1084  G4ParticleDefinition * Quark;
1085 
1086  if(string->GetLeftParton()->GetParticleSubType()== "quark")
1087  {
1088  Quark =string->GetLeftParton();
1089  Di_Quark=string->GetRightParton();
1090  } else
1091  {
1092  Quark =string->GetRightParton();
1093  Di_Quark=string->GetLeftParton();
1094  }
1095 
1096  G4int IDquark =Quark->GetPDGEncoding();
1097  G4int AbsIDquark =std::abs(IDquark);
1098  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1099  G4int AbsIDdi_quark=std::abs(IDdi_quark);
1100  G4int Di_q1=AbsIDdi_quark/1000;
1101  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1102  //G4cout<<"IDs "<<IDdi_quark<<" "<<IDquark<<G4endl;
1103 
1104  G4int SignDiQ= 1;
1105  if(IDdi_quark < 0) SignDiQ=-1;
1106 
1107  NumberOf_FS=0;
1108  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1109  {
1110  G4int SignQ;
1111  if(IDquark > 0)
1112  { SignQ=-1;
1113  if(IDquark == 2) SignQ= 1;
1114  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1115  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1116  } else
1117  {
1118  SignQ= 1;
1119  if(IDquark == -2) SignQ=-1;
1120  if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1121  if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1122  }
1123 
1124  if(AbsIDquark == ProdQ) SignQ= 1;
1125 
1126  //G4cout<<G4endl;
1127  //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl;
1128 
1129  G4int StateQ=0;
1130  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1131  {
1132  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1133  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1134  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1135  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1136  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1137 
1138  //G4cout<<"Meson "<<LeftHadron->GetParticleName()<<G4endl;
1139 
1140  G4int StateDiQ=0;
1141  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1142  {
1143  //G4cout<<G4endl<<"Di_q1 Di_q2 ProdQ StateDiQ "<<BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1144  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1145  Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1146  G4double RightHadronMass=RightHadron->GetPDGMass();
1147 
1148  //G4cout<<"Baryon "<<RightHadron->GetParticleName()<<G4endl;
1149 
1150  //G4cout<<"StringMass LeftHadronMass RightHadronMass "<<StringMass<<" "<<LeftHadronMass<<" "<< RightHadronMass<<G4endl;
1151 
1152  if(StringMass >= LeftHadronMass + RightHadronMass)
1153  {
1154  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1155  sqr(RightHadronMass));
1156  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1157  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1158  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1159  Prob_QQbar[ProdQ-1];
1160 
1161  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1162  FS_RightHadron[NumberOf_FS]= RightHadron;
1163 
1164  //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl;
1165  //G4cout<<"++++++++++++++++++++++++++++++++"<<G4endl;
1166  NumberOf_FS++;
1167 
1168  if(NumberOf_FS > 34)
1169  {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;}
1170  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1171 
1172  StateDiQ++;
1173  //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1174  } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0);
1175 
1176  StateQ++;
1177  } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
1178  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1179 
1180  return true;
1181 }
1182 
1183 //----------------------------------------------------------------------------------------
1184 G4bool G4LundStringFragmentation::
1185 Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
1186  G4ParticleDefinition * & LeftHadron,
1187  G4ParticleDefinition * & RightHadron)
1188 {
1189  G4double StringMass = string->Mass();
1190  G4double StringMassSqr= sqr(StringMass);
1191 
1192  G4ParticleDefinition * Quark;
1193  G4ParticleDefinition * Anti_Quark;
1194 
1195  if(string->GetLeftParton()->GetPDGEncoding()>0)
1196  {
1197  Quark =string->GetLeftParton();
1198  Anti_Quark=string->GetRightParton();
1199  } else
1200  {
1201  Quark =string->GetRightParton();
1202  Anti_Quark=string->GetLeftParton();
1203  }
1204 
1205  //G4cout<<" Quarks "<<Quark->GetParticleName()<<" "<<Anti_Quark->GetParticleName()<<G4endl;
1206  G4int IDquark =Quark->GetPDGEncoding();
1207  G4int AbsIDquark =std::abs(IDquark);
1208  G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1209  G4int AbsIDanti_quark=std::abs(IDanti_quark);
1210 
1211  NumberOf_FS=0;
1212  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1213  {
1214  //G4cout<<"ProdQ "<<ProdQ<<G4endl;
1215  G4int SignQ=-1;
1216  if(IDquark == 2) SignQ= 1;
1217  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1218  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1219  if(IDquark == ProdQ) SignQ= 1;
1220 
1221  G4int SignAQ= 1;
1222  if(IDanti_quark == -2) SignAQ=-1;
1223  if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
1224  if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
1225  if(AbsIDanti_quark == ProdQ) SignAQ= 1;
1226 
1227  G4int StateQ=0;
1228  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1229  {
1230  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1231  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1232  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1233  StateQ++;
1234 
1235  G4int StateAQ=0;
1236  do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0);
1237  {
1238  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1239  Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1240  G4double RightHadronMass=RightHadron->GetPDGMass();
1241  StateAQ++;
1242 
1243  if(StringMass >= LeftHadronMass + RightHadronMass)
1244  {
1245  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1246  sqr(RightHadronMass));
1247  //FS_Psqr=1.;
1248  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1249  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1250  MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1251  Prob_QQbar[ProdQ-1];
1252 
1253  if(string->GetLeftParton()->GetPDGEncoding()>0)
1254  {
1255  FS_LeftHadron[NumberOf_FS] = RightHadron;
1256  FS_RightHadron[NumberOf_FS]= LeftHadron;
1257  } else
1258  {
1259  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1260  FS_RightHadron[NumberOf_FS]= RightHadron;
1261  }
1262  NumberOf_FS++;
1263  //G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<" ";
1264  //G4cout<<"Masses "<<StringMass<<" "<<LeftHadronMass<<" "<<RightHadronMass<<" "<<NumberOf_FS-1<<G4endl; //FS_Psqr<<G4endl;
1265 
1266  if(NumberOf_FS > 34)
1267  {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;}
1268  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1269  } while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0);
1270  } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
1271  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1272 
1273  return true;
1274 }
1275 
1276 //----------------------------------------------------------------------------------------------------------
1277 G4int G4LundStringFragmentation::SampleState(void)
1278 {
1279  G4double SumWeights=0.;
1280  for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1281 
1282  G4double ksi=G4UniformRand();
1283  G4double Sum=0.;
1284  G4int indexPosition = 0;
1285 
1286  for(G4int i=0; i<NumberOf_FS; i++)
1287  {
1288  Sum+=(FS_Weight[i]/SumWeights);
1289  indexPosition=i;
1290  if(Sum >= ksi) break;
1291  }
1292  return indexPosition;
1293 }
1294 
G4LundStringFragmentation::FragmentString
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
Definition: G4LundStringFragmentation.cc:492
G4FragmentingString::GetRightParton
G4ParticleDefinition * GetRightParton(void) const
Definition: G4FragmentingString.hh:143
G4FragmentingString.hh
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:620
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:237
G4VLongitudinalStringDecay::CPExcited
G4ExcitedString * CPExcited(const G4ExcitedString &string)
Definition: G4VLongitudinalStringDecay.cc:314
Pos
ush Pos
Definition: deflate.h:89
G4KineticTrack::SetFormationTime
void SetFormationTime(G4double aFormationTime)
Definition: G4KineticTrack.hh:230
G4VLongitudinalStringDecay::SetStrangenessSuppression
void SetStrangenessSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:525
G4HadronBuilder::Build
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition: G4HadronBuilder.cc:53
G4FragmentingString::Mass
G4double Mass() const
Definition: G4FragmentingString.cc:265
CLHEP::Hep3Vector::getY
double getY() const
G4VLongitudinalStringDecay::SampleQuarkPt
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
Definition: G4VLongitudinalStringDecay.cc:468
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:145
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::setZ
void setZ(double)
G4FragmentingString::GetDecayParton
G4ParticleDefinition * GetDecayParton() const
Definition: G4FragmentingString.hh:131
G4VLongitudinalStringDecay::SetStringTensionParameter
void SetStringTensionParameter(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:628
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:137
G4VLongitudinalStringDecay::SetDiquarkBreakProbability
void SetDiquarkBreakProbability(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:547
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4VLongitudinalStringDecay::LightFragmentationTest
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
Definition: G4VLongitudinalStringDecay.cc:151
G4ExcitedString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4ExcitedString.hh:161
G4KineticTrack::GetFormationTime
G4double GetFormationTime() const
Definition: G4KineticTrack.hh:225
CLHEP::HepLorentzVector::mag
double mag() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4ExcitedString::TransformToAlignedCms
G4LorentzRotation TransformToAlignedCms()
Definition: G4ExcitedString.hh:215
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
G4VLongitudinalStringDecay::FindParticle
G4ParticleDefinition * FindParticle(G4int Encoding)
Definition: G4VLongitudinalStringDecay.cc:281
G4KineticTrack::SetPosition
void SetPosition(const G4ThreeVector aPosition)
Definition: G4KineticTrack.hh:242
Randomize.hh
G4FragmentingString::GetDecayDirection
G4int GetDecayDirection() const
Definition: G4FragmentingString.cc:184
G4KineticTrack::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4KineticTrack.hh:258
G4Quarks.hh
G4LundStringFragmentation.hh
G4FragmentingString::FourQuarkString
G4bool FourQuarkString(void) const
Definition: G4FragmentingString.cc:194
G4VLongitudinalStringDecay::scalarMesonMix
std::vector< G4double > scalarMesonMix
Definition: G4VLongitudinalStringDecay.hh:190
G4PhysicalConstants.hh
G4DiQuarks.hh
C1
#define C1
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:95
CLHEP::HepLorentzVector::mt2
double mt2() const
G4ExcitedString::GetTimeOfCreation
G4double GetTimeOfCreation() const
Definition: G4ExcitedString.hh:137
CLHEP::Hep3Vector::mag2
double mag2() const
test.z
tuple z
Definition: test.py:28
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
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:248
G4VLongitudinalStringDecay::CalculateHadronTimePosition
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
Definition: G4VLongitudinalStringDecay.cc:485
python.hepunit.c_light
float c_light
Definition: hepunit.py:257
G4VLongitudinalStringDecay::Splitup
G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
Definition: G4VLongitudinalStringDecay.cc:323
G4VLongitudinalStringDecay::SetDiquarkSuppression
void SetDiquarkSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:536
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
Generated on Sat Dec 14 2013 14:34:59 for Geant4_10 by   doxygen 1.8.5