source: trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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