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