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

Last change on this file since 1344 was 1340, checked in by garnier, 15 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.