Geant4  10.00.p03
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 81884 2014-06-06 13:03:02Z 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 // -----------------------------------------------------------------------------
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 
47 {
48 // ------ For estimation of a minimal string mass ---------------
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);
63  SetDiquarkSuppression(0.087); // Uzhi 18.05.2012
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
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
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
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 // --------------------------------------------------------------
403 {}
404 
405 
406 //--------------------------------------------------------------------------------------
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 //--------------------------------------------------------------------------------------
486  const G4double aValue)
487 {
488  MinimalStringMass2=aValue * aValue;
489 }
490 
491 //--------------------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------
583 {
584  SetMinimalStringMass(string);
585  // return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
586  return MinimalStringMass < string->Get4Momentum().mag(); // 21.07.2010
587 }
588 
589 //----------------------------------------------------------------------------------------
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() -
600  }
601 }
602 
603 //----------------------------------------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------------------------------
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 //-----------------------------------------------------------------------------
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 //-----------------------------------------------------------------------------------------
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 //------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------------
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;
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;
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  //M. Novak: = was removed.
1043  // 1. it can cause numerical problems
1044  // 2. weight of the corresponding FS is zero (in theory)
1045  if(StringMass > LeftHadronMass + RightHadronMass)
1046  {
1047  if ( NumberOf_FS > 34 ) {
1049  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1050  G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
1051  "HAD_LUND_001", JustWarning, ed );
1052  NumberOf_FS = 34;
1053  }
1054 
1055  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1056  sqr(RightHadronMass));
1057  //FS_Psqr=1.;
1058  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
1059  BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
1060  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1061  Prob_QQbar[ProdQ-1];
1062 
1063  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1064  FS_RightHadron[NumberOf_FS]= RightHadron;
1065 
1066  //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;
1067  //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl;
1068  NumberOf_FS++;
1069  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1070 
1071  StateDiQ++;
1072  //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1073  } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0);
1074 
1075  StateADiQ++;
1076  } while(Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0);
1077  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1078 
1079  return true;
1080 }
1081 
1082 //----------------------------------------------------------------------------------------
1085  G4ParticleDefinition * & LeftHadron,
1086  G4ParticleDefinition * & RightHadron)
1087 {
1088  G4double StringMass = string->Mass();
1089  G4double StringMassSqr= sqr(StringMass);
1090 
1091  G4ParticleDefinition * Di_Quark;
1092  G4ParticleDefinition * Quark;
1093 
1094  if(string->GetLeftParton()->GetParticleSubType()== "quark")
1095  {
1096  Quark =string->GetLeftParton();
1097  Di_Quark=string->GetRightParton();
1098  } else
1099  {
1100  Quark =string->GetRightParton();
1101  Di_Quark=string->GetLeftParton();
1102  }
1103 
1104  G4int IDquark =Quark->GetPDGEncoding();
1105  G4int AbsIDquark =std::abs(IDquark);
1106  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1107  G4int AbsIDdi_quark=std::abs(IDdi_quark);
1108  G4int Di_q1=AbsIDdi_quark/1000;
1109  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1110  //G4cout<<"IDs "<<IDdi_quark<<" "<<IDquark<<G4endl;
1111 
1112  G4int SignDiQ= 1;
1113  if(IDdi_quark < 0) SignDiQ=-1;
1114 
1115  NumberOf_FS=0;
1116  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1117  {
1118  G4int SignQ;
1119  if(IDquark > 0)
1120  { SignQ=-1;
1121  if(IDquark == 2) SignQ= 1;
1122  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1123  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1124  } else
1125  {
1126  SignQ= 1;
1127  if(IDquark == -2) SignQ=-1;
1128  if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1129  if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1130  }
1131 
1132  if(AbsIDquark == ProdQ) SignQ= 1;
1133 
1134  //G4cout<<G4endl;
1135  //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl;
1136 
1137  G4int StateQ=0;
1138  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1139  {
1140  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1141  //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1142  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1143  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1144  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1145 
1146  //G4cout<<"Meson "<<LeftHadron->GetParticleName()<<G4endl;
1147 
1148  G4int StateDiQ=0;
1149  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1150  {
1151  //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;
1152  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1153  Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1154  G4double RightHadronMass=RightHadron->GetPDGMass();
1155 
1156  //G4cout<<"Baryon "<<RightHadron->GetParticleName()<<G4endl;
1157 
1158  //G4cout<<"StringMass LeftHadronMass RightHadronMass "<<StringMass<<" "<<LeftHadronMass<<" "<< RightHadronMass<<G4endl;
1159 
1160  //M. Novak: = was removed.
1161  // 1. it can cause numerical problems
1162  // 2. weight of the corresponding FS is zero (in theory)
1163  if(StringMass > LeftHadronMass + RightHadronMass)
1164  {
1165  if ( NumberOf_FS > 34 ) {
1167  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1168  G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1169  "HAD_LUND_002", JustWarning, ed );
1170  NumberOf_FS = 34;
1171  }
1172 
1173  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1174  sqr(RightHadronMass));
1175  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1176  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1177  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1178  Prob_QQbar[ProdQ-1];
1179 
1180  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1181  FS_RightHadron[NumberOf_FS]= RightHadron;
1182 
1183  //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl;
1184  //G4cout<<"++++++++++++++++++++++++++++++++"<<G4endl;
1185  NumberOf_FS++;
1186  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1187 
1188  StateDiQ++;
1189  //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
1190  } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0);
1191 
1192  StateQ++;
1193  } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
1194  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1195 
1196  return true;
1197 }
1198 
1199 //----------------------------------------------------------------------------------------
1202  G4ParticleDefinition * & LeftHadron,
1203  G4ParticleDefinition * & RightHadron)
1204 {
1205  G4double StringMass = string->Mass();
1206  G4double StringMassSqr= sqr(StringMass);
1207 
1208  G4ParticleDefinition * Quark;
1209  G4ParticleDefinition * Anti_Quark;
1210 
1211  if(string->GetLeftParton()->GetPDGEncoding()>0)
1212  {
1213  Quark =string->GetLeftParton();
1214  Anti_Quark=string->GetRightParton();
1215  } else
1216  {
1217  Quark =string->GetRightParton();
1218  Anti_Quark=string->GetLeftParton();
1219  }
1220 
1221  //G4cout<<" Quarks "<<Quark->GetParticleName()<<" "<<Anti_Quark->GetParticleName()<<G4endl;
1222  G4int IDquark =Quark->GetPDGEncoding();
1223  G4int AbsIDquark =std::abs(IDquark);
1224  G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1225  G4int AbsIDanti_quark=std::abs(IDanti_quark);
1226 
1227  NumberOf_FS=0;
1228  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1229  {
1230  //G4cout<<"ProdQ "<<ProdQ<<G4endl;
1231  G4int SignQ=-1;
1232  if(IDquark == 2) SignQ= 1;
1233  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1234  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1235  if(IDquark == ProdQ) SignQ= 1;
1236 
1237  G4int SignAQ= 1;
1238  if(IDanti_quark == -2) SignAQ=-1;
1239  if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
1240  if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
1241  if(AbsIDanti_quark == ProdQ) SignAQ= 1;
1242 
1243  G4int StateQ=0;
1244  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1245  {
1246  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1247  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1248  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1249  StateQ++;
1250 
1251  G4int StateAQ=0;
1252  do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0);
1253  {
1254  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1255  Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1256  G4double RightHadronMass=RightHadron->GetPDGMass();
1257  StateAQ++;
1258 
1259  //M. Novak: = was removed.
1260  // 1. it can cause numerical problems
1261  // 2. weight of the corresponding FS is zero (in theory)
1262  if(StringMass > LeftHadronMass + RightHadronMass)
1263  {
1264  if ( NumberOf_FS > 34 ) {
1266  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1267  G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1268  "HAD_LUND_003", JustWarning, ed );
1269  NumberOf_FS = 34;
1270  }
1271 
1272  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1273  sqr(RightHadronMass));
1274  //FS_Psqr=1.;
1275  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1276  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1277  MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1278  Prob_QQbar[ProdQ-1];
1279 
1280  if(string->GetLeftParton()->GetPDGEncoding()>0)
1281  {
1282  FS_LeftHadron[NumberOf_FS] = RightHadron;
1283  FS_RightHadron[NumberOf_FS]= LeftHadron;
1284  } else
1285  {
1286  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1287  FS_RightHadron[NumberOf_FS]= RightHadron;
1288  }
1289  NumberOf_FS++;
1290  //G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<" ";
1291  //G4cout<<"Masses "<<StringMass<<" "<<LeftHadronMass<<" "<<RightHadronMass<<" "<<NumberOf_FS-1<<G4endl; //FS_Psqr<<G4endl;
1292  } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
1293  } while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0);
1294  } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
1295  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1296 
1297  return true;
1298 }
1299 
1300 //----------------------------------------------------------------------------------------------------------
1302 {
1303  if ( NumberOf_FS > 34 ) {
1305  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1306  G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1307  NumberOf_FS = 34;
1308  }
1309 
1310  G4double SumWeights=0.;
1311  for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1312 
1313  G4double ksi=G4UniformRand();
1314  G4double Sum=0.;
1315  G4int indexPosition = 0;
1316 
1317  for(G4int i=0; i<NumberOf_FS; i++)
1318  {
1319  Sum+=(FS_Weight[i]/SumWeights);
1320  indexPosition=i;
1321  if(Sum >= ksi) break;
1322  }
1323  return indexPosition;
1324 }
1325 
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
G4ParticleDefinition * GetRightParton(void) const
static const double MeV
Definition: G4SIunits.hh:193
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4bool Quark_AntiQuark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
CLHEP::HepLorentzRotation G4LorentzRotation
const G4ThreeVector & GetPosition() const
G4ExcitedString * CPExcited(const G4ExcitedString &string)
G4double z
Definition: TRTMaterials.hh:39
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
const G4double pi
void SetFormationTime(G4double aFormationTime)
void SetStrangenessSuppression(G4double aValue)
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
G4bool Quark_Diquark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
G4double Mass() const
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
std::vector< G4double > vectorMesonMix
const G4String & GetParticleSubType() const
virtual void SetMassCut(G4double aValue)
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition * GetDecayParton() const
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
void SetStringTensionParameter(G4double aValue)
const G4ThreeVector & GetPosition() const
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
G4ParticleDefinition * GetLeftParton(void) const
void SetDiquarkBreakProbability(G4double aValue)
void SetMinimalStringMass2(const G4double aValue)
#define G4UniformRand()
Definition: Randomize.hh:87
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
G4LorentzVector Get4Momentum() const
G4bool Loop_toFragmentString(G4ExcitedString *&theStringInCMS, G4KineticTrackVector *&LeftVector, G4KineticTrackVector *&RightVector)
G4double GetFormationTime() const
void SetMinimalStringMass(const G4FragmentingString *const string)
bool G4bool
Definition: G4Types.hh:79
G4LorentzRotation TransformToAlignedCms()
G4ParticleDefinition * FindParticle(G4int Encoding)
void SetPosition(const G4ThreeVector aPosition)
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
static const double GeV
Definition: G4SIunits.hh:196
G4int GetDecayDirection() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4double lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
G4bool FourQuarkString(void) const
std::vector< G4double > scalarMesonMix
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define C1
G4double GetPDGMass() const
G4bool Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
static G4ParticleTable * GetParticleTable()
G4double GetTimeOfCreation() const
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * FS_LeftHadron[35]
G4bool Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
T sqr(const T &x)
Definition: templates.hh:145
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
static const double fermi
Definition: G4SIunits.hh:93
G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
G4ParticleDefinition * FS_RightHadron[35]
void SetDiquarkSuppression(G4double aValue)
CLHEP::HepLorentzVector G4LorentzVector