| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // History: first implementation, Maxim K
32 // -------------------------------------------
33 #include "G4LundStringFragmentation.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 #include "G4FragmentingString.hh"
38 #include "G4DiQuarks.hh"
39 #include "G4Quarks.hh"
40 #include "G4HadronicParameters.hh"
41 #include "G4Exp.hh"
42 #include "G4Pow.hh"
43
44 //#define debug_LUNDfragmentation
45
46 // Class G4LundStringFragmentation
47 //********************************************
48
49 G4LundStringFragmentation::G4LundStringFragmen
50 : G4VLongitudinalStringDecay("LundStringFrag
51 {
52 SetMassCut(210.*MeV); // Mpi + Delta
53 // For ProduceOneH
54 // that no one pi-
55 SigmaQT = 0.435 * GeV;
56 Tmt = 190.0 * MeV;
57
58 SetStringTensionParameter(1.*GeV/fermi);
59 SetDiquarkBreakProbability(0.3);
60
61 SetStrangenessSuppression((1.0 - 0.12)/2.0
62 SetDiquarkSuppression(0.07);
63
64 // Check if charmed and bottom hadrons are
65 // set the non-zero probabilities for c-cb
66 // else set them to 0.0. If these probabil
67 // hadrons can't/can be created during the
68 // (i.e. not heavy) projectile hadron nucl
69 if ( G4HadronicParameters::Instance()->Ena
70 SetProbCCbar(0.0002); // According to O
71 SetProbBBbar(5.0e-5); // According to O
72 } else {
73 SetProbCCbar(0.0);
74 SetProbBBbar(0.0);
75 }
76
77 SetMinMasses(); // For treating of small
78 }
79
80 //--------------------------------------------
81
82 G4KineticTrackVector* G4LundStringFragmentatio
83 {
84 // Can no longer modify Parameters for Fragm
85
86 PastInitPhase=true;
87
88 G4FragmentingString aString(theString);
89 SetMinimalStringMass(&aString);
90
91 #ifdef debug_LUNDfragmentation
92 G4cout<<G4endl<<"LUND StringFragmentat
93 G4cout<<G4endl<<"LUND StringFragm: Str
94 <<theString.Get4Momentum
95 <<"4Mom "<<theString.Get
96 <<"---------------------
97 G4cout<<"String ends Direct "<<theStri
98 <<theStri
99 <<theStri
100 G4cout<<"Left mom "<<theString.GetLef
101 G4cout<<"Right mom "<<theString.GetRig
102 G4cout<<"Check for Fragmentation "<<G4
103 #endif
104
105 G4KineticTrackVector * LeftVector(0);
106
107 if (!aString.IsAFourQuarkString() && !IsItFr
108 {
109 #ifdef debug_LUNDfragmentation
110 G4cout<<"Non fragmentable - th
111 #endif
112 // SetMassCut(210.*MeV); // F
113 // t
114
115 G4double Mcut = GetMassCut();
116 SetMassCut(10000.*MeV);
117 LeftVector=ProduceOneHadron(&theString);
118 SetMassCut(Mcut);
119
120 if ( LeftVector )
121 {
122 if ( LeftVector->size() > 0)
123 {
124 LeftVector->operator[](0)->SetForm
125 LeftVector->operator[](0)->SetPosi
126 }
127 if (LeftVector->size() > 1)
128 {
129 // 2 hadrons created from qq-qqbar
130 LeftVector->operator[](1)->SetFormationT
131 LeftVector->operator[](1)->SetPosition(t
132 }
133 }
134 return LeftVector;
135 }
136
137 #ifdef debug_LUNDfragmentation
138 G4cout<<"The string will be fragmented
139 #endif
140
141 // The string can fragment. At least two par
142 LeftVector =new G4KineticTrackVec
143 G4KineticTrackVector * RightVector=new G4Kin
144
145 G4bool success = Loop_toFragmentString
146
147 if ( ! success )
148 {
149 std::for_each(LeftVector->begin(), LeftVec
150 LeftVector->clear();
151 std::for_each(RightVector->begin(), RightV
152 delete RightVector;
153 return LeftVector;
154 }
155
156 // Join Left- and RightVector into LeftVecto
157 while (!RightVector->empty())
158 {
159 LeftVector->push_back(RightVector->back())
160 RightVector->erase(RightVector->end()-1);
161 }
162 delete RightVector;
163
164 return LeftVector;
165 }
166
167 //--------------------------------------------
168
169 G4bool G4LundStringFragmentation::IsItFragment
170 {
171 SetMinimalStringMass(string);
172 //G4cout<<"MinM StrM "<<MinimalStringM
173
174 return std::abs(MinimalStringMass) < string-
175
176 //MinimalStringMass is negative and la
177 }
178
179 //--------------------------------------------
180
181 G4bool G4LundStringFragmentation::Loop_toFragm
182
183
184 {
185 #ifdef debug_LUNDfragmentation
186 G4cout<<"Loop_toFrag "<<theString.GetL
187 <<theString.GetL
188 <<" "<<theString.GetR
189 <<theString.GetR
190 <<"Direction "<<theString.GetD
191 #endif
192
193 G4LorentzRotation toCmsI, toObserverFr
194
195 G4bool final_success=false;
196 G4bool inner_success=true;
197
198 G4int attempt=0;
199
200 while ( ! final_success && attempt++ < Strin
201 { // If the string fragmentation does
202 // repeat the fragmentation.
203
204 G4FragmentingString *currentSt
205 toCmsI = currentString->Transf
206 toObserverFrameI = toCmsI.inve
207
208 G4LorentzRotation toCms, toObserverFrame;
209
210 //G4cout<<"Main loop start whilecounter "<
211
212 // Cleaning up the previously produced had
213 std::for_each(LeftVector->begin(), LeftVec
214 LeftVector->clear();
215 std::for_each(RightVector->begin(), RightV
216 RightVector->clear();
217
218 // Main fragmentation loop until the strin
219 inner_success=true; // set false on failu
220 const G4int maxNumberOfLoops =
221 G4int loopCounter = -1;
222
223 while ( (! StopFragmenting(currentString))
224 { // Split current string into hadro
225 #ifdef debug_LUNDfragm
226 G4cout<<"The string wi
227 //G4cout<<"1 "<<curren
228 #endif
229 G4FragmentingString *newString=0; // us
230
231 toCms=currentString->TransformToAlignedC
232 toObserverFrame= toCms
233
234 #ifdef debug_LUNDfragm
235 //G4cout<<"CMS Left m
236 //G4cout<<"CMS Right m
237 //G4cout<<"CMS String
238 #endif
239
240 G4KineticTrack * Hadron=Splitup(currentS
241
242 if ( Hadron != 0 ) // Store the hadron
243 {
244 #ifdef debug_L
245 G4cout<<"Hadro
246 //G4cout<<"2 "
247 #endif
248
249 Hadron->Set4Momentum(toObserverFrame*H
250
251 G4double TimeOftheStringCreation=theSt
252 G4ThreeVector PositionOftheStringCreat
253
254 G4LorentzVector Coordinate(Hadron->Get
255 G4LorentzVector Momentum = toObserverF
256 Hadron->SetFormationTime(TimeOftheStri
257 G4ThreeVector aPosition(Momentum.vect(
258 Hadron->SetPosition(PositionOftheStrin
259
260 // Open to pro
261 if ( currentString->GetDecayDirection(
262 {
263 LeftVector->push_back(Hadron);
264 } else
265 {
266 RightVector->push_back(Hadron);
267 }
268 delete currentString;
269 currentString=newString;
270 } else {
271 if ( newString ) del
272 }
273
274 currentString->Lorentz
275 };
276
277 if ( loopCounter >= maxNumberO
278 inner_success=false;
279 }
280
281 // Split remaining string into 2 final had
282 #ifdef debug_LUNDfragmentation
283 if (inner_success) G4cout<<"Sp
284 #endif
285
286 if ( inner_success && SplitLast(currentStr
287 {
288 final_success = true;
289 }
290
291 delete currentString;
292 } // End of the loop where we try to fragme
293
294 G4int sign = +1;
295 if ( theString.GetDirection() < 0 ) si
296 for ( unsigned int hadronI = 0; hadron
297 G4LorentzVector Tmp = LeftVector->o
298 Tmp.setZ(sign*Tmp.getZ());
299 Tmp *= toObserverFrameI;
300 LeftVector->operator[](hadronI)->Se
301 }
302 for ( unsigned int hadronI = 0; hadron
303 G4LorentzVector Tmp = RightVector->
304 Tmp.setZ(sign*Tmp.getZ());
305 Tmp *= toObserverFrameI;
306 RightVector->operator[](hadronI)->S
307 }
308
309 return final_success;
310 }
311
312 //--------------------------------------------
313
314 G4bool G4LundStringFragmentation::StopFragment
315 {
316 SetMinimalStringMass(string);
317
318 if ( MinimalStringMass < 0.) return true;
319
320 if (string->IsAFourQuarkString())
321 {
322 return G4UniformRand() < G4Exp(-0.0005*(st
323 } else {
324
325 if (MinimalStringMass < 0.0 )
326
327 G4bool Result = G4UniformRand() <
328 G4Exp(-0.66e-6*(string->Mass()*string-
329 // G4bool Result = string->Mas
330
331 #ifdef debug_LUNDfragmentation
332 G4cout<<"StopFragmenting Minim
333 <<" "<<string->Mass()<<G
334 G4cout<<"StopFragmenting - Yes
335 #endif
336 return Result;
337 }
338 }
339
340 //--------------------------------------------
341
342 G4KineticTrack * G4LundStringFragmentation::Sp
343 G4Fragmentin
344 {
345 #ifdef debug_LUNDfragmentation
346 G4cout<<G4endl;
347 G4cout<<"Start SplitUP ================
348 G4cout<<"String partons: " <<string->Ge
349 <<string->Ge
350 <<"Direction " <<string->Ge
351 #endif
352
353 //... random choice of string end to us
354 G4int SideOfDecay = (G4UniformRand() <
355 if (SideOfDecay < 0)
356 {
357 string->SetLeftPartonStable();
358 } else
359 {
360 string->SetRightPartonStable();
361 }
362
363 G4ParticleDefinition *newStringEnd;
364 G4ParticleDefinition * HadronDefinition
365
366 G4double StringMass=string->Mass();
367
368 G4double ProbDqADq = GetDiquarkSuppress
369 G4double ProbSaS = 1.0 - 2.0 * GetStr
370
371 #ifdef debug_LUNDfragmentation
372 G4cout<<"StrMass DiquarkSuppression
373 #endif
374
375 G4int NumberOfpossibleBaryons = 2;
376
377 if (string->GetLeftParton()->GetParticl
378 if (string->GetRightParton()->GetPartic
379
380 G4double ActualProb = ProbDqADq ;
381 ActualProb *= (1.0-G4Pow::GetInstance()
382 if(ActualProb <0.0) ActualProb = 0.;
383
384 SetDiquarkSuppression(ActualProb);
385
386 G4double Mth = 1250.0;
387 if ( NumberOfpossibleBaryons == 3 ){Mth
388 else if ( NumberOfpossibleBaryons == 4
389 else {}
390
391 ActualProb = ProbSaS;
392 ActualProb *= (1.0 - G4Pow::GetInstance
393 if ( ActualProb < 0.0 ) ActualProb = 0.
394 SetStrangenessSuppression((1.0-ActualPr
395
396 #ifdef debug_LUNDfragmentation
397 G4cout<<"StrMass DiquarkSuppression cor
398 #endif
399
400 if (string->DecayIsQuark())
401 {
402 HadronDefinition= QuarkSplitup(strin
403 } else {
404 HadronDefinition= DiQuarkSplitup(str
405 }
406
407 SetDiquarkSuppression(ProbDqADq);
408 SetStrangenessSuppression((1.0-ProbSaS)
409
410 if ( HadronDefinition == NULL ) { G4Kin
411
412 #ifdef debug_LUNDfragmentation
413 G4cout<<"The parton "<<string->GetDecay
414 <<" produces hadron "<<HadronDefi
415 <<" and is transformed to "<<newS
416 G4cout<<"The side of the string decay L
417 #endif
418 // create new String from old, ie. keep
419
420 if ( newString ) delete newString;
421
422 newString=new G4FragmentingString(*stri
423
424 #ifdef debug_LUNDfragmentation
425 G4cout<<"An attempt to determine its en
426 #endif
427 G4LorentzVector* HadronMomentum=SplitEa
428
429 delete newString; newString=0;
430
431 G4KineticTrack * Hadron =0;
432 if ( HadronMomentum != 0 ) {
433
434 #ifdef debug_LUNDfragmentation
435 G4cout<<"The attempt was successful
436 #endif
437 G4ThreeVector Pos;
438 Hadron = new G4KineticTrack(HadronDefinit
439
440 if ( newString ) delete newString;
441
442 newString=new G4FragmentingString(*string
443 HadronMomentum);
444 delete HadronMomentum;
445 }
446 else
447 {
448 #ifdef debug_LUNDfragmentation
449 G4cout<<"The attempt was not succes
450 #endif
451 }
452
453 #ifdef debug_LUNDfragmentation
454 G4cout<<"End SplitUP (G4VLongitudinalSt
455 #endif
456
457 return Hadron;
458 }
459
460 //--------------------------------------------
461
462 G4ParticleDefinition * G4LundStringFragmentati
463
464 {
465 G4double StrSup=GetStrangeSuppress();
466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2
467
468 //... can Diquark break or not?
469 if (G4UniformRand() < DiquarkBreakProb ){
470
471 //... Diquark break
472 G4int stableQuarkEncoding = decay->GetPD
473 G4int decayQuarkEncoding = (decay->GetPD
474 if (G4UniformRand() < 0.5)
475 {
476 G4int Swap = stableQuarkEncoding;
477 stableQuarkEncoding = decayQuarkEncod
478 decayQuarkEncoding = Swap;
479 }
480
481 G4int IsParticle=(decayQuarkEncoding>0)
482
483 SetStrangenessSuppression((1.0-ProbQQbar
484 pDefPair QuarkPair = CreatePartonPair(Is
485 SetStrangenessSuppression((1.0-StrSup)/2
486
487 //... Build new Diquark
488 G4int QuarkEncoding=QuarkPair.second->Ge
489 G4int i10 = std::max(std::abs(QuarkEnco
490 G4int i20 = std::min(std::abs(QuarkEnco
491 G4int spin = (i10 != i20 && G4UniformRan
492 G4int NewDecayEncoding = -1*IsParticle*(
493 created = FindParticle(NewDecayEncoding)
494 G4ParticleDefinition * decayQuark=FindPa
495 G4ParticleDefinition * had=hadronizer->B
496 StrangeSuppress=StrSup;
497
498 return had;
499
500 } else {
501 //... Diquark does not break
502
503 G4int IsParticle=(decay->GetPDGEncoding(
504
505 StrangeSuppress=(1.0 - ProbQQbar)/2.0;
506 pDefPair QuarkPair = CreatePartonPair(Is
507
508 created = QuarkPair.second;
509
510 G4ParticleDefinition * had=hadronizer->B
511 StrangeSuppress=StrSup;
512
513 return had;
514 }
515 }
516
517 //--------------------------------------------
518
519 G4LorentzVector * G4LundStringFragmentation::S
520
521
522 {
523 G4LorentzVector String4Momentum=string->Get4
524 G4double StringMT2=string->MassT2();
525 G4double StringMT =std::sqrt(StringMT2);
526
527 G4double HadronMass = pHadron->GetPDGMass();
528 SetMinimalStringMass(newString);
529
530 if ( MinimalStringMass < 0.0 ) return n
531
532 #ifdef debug_LUNDfragmentation
533 G4cout<<G4endl<<"Start LUND SplitEandP
534 G4cout<<"String 4 mom, String M and Mt
535 <<" "<<std::sqrt(StringMT2)<<G4e
536 G4cout<<"Hadron "<<pHadron->GetParticl
537 G4cout<<"HadM MinimalStringMassLeft St
538 <<String4Momentum.mag()<<" "<<Ha
539 #endif
540
541 if ((HadronMass + MinimalStringMass > string
542 {
543 #ifdef debug_LUNDfragmentation
544 G4cout<<"Mass of the string is not s
545 #endif
546 return 0;
547 } // have to start all over!
548
549 String4Momentum.setPz(0.);
550 G4ThreeVector StringPt=String4Momentum.vect(
551 StringPt.setZ(0.);
552
553 // calculate and assign hadron transverse mo
554 G4ThreeVector HadronPt , RemSysPt;
555 G4double HadronMassT2, ResidualMassT2;
556 G4double HadronMt, Pt, Pt2, phi;
557
558 G4double TmtCur = Tmt;
559
560 if ( (string->GetDecayParton()->GetPar
561 (pHadron->GetBaryonNumber() != 0)
562 TmtCur = Tmt*0.37; // q
563 } else if ( (string->GetDecayParton()-
564 (pHadron->GetBaryonNumber(
565 //TmtCur = Tmt;
566 } else if ( (string->GetDecayParton()->GetPa
567 (pHadron->GetBaryonNumber(
568 //TmtCur = Tmt*0.89;
569 } else if ( (string->GetDecayParton()-
570 (pHadron->GetBaryonNumber(
571 TmtCur = Tmt*1.35;
572 }
573
574 //... sample Pt of the hadron
575 G4int attempt=0;
576 do
577 {
578 attempt++; if (attempt > StringLoopI
579
580 HadronMt = HadronMass - TmtCur*G4Log
581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st
582 phi = 2.*pi*G4UniformRand();
583 HadronPt = G4ThreeVector( Pt*std::co
584 RemSysPt = StringPt - HadronPt;
585 HadronMassT2 = sqr(HadronMass) + Had
586 ResidualMassT2=sqr(MinimalStringMass
587
588 } while (std::sqrt(HadronMassT2) + std
589
590 //... sample z to define hadron longitudina
591 //... but first check the available phase sp
592
593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2
594 4*HadronMassT2 * ResidualMassT2)/4./Stri
595
596 if (Pz2 < 0 ) {return 0;} // have t
597
598 //... then compute allowed z region z_min <
599
600 G4double Pz = std::sqrt(Pz2);
601 G4double zMin = (std::sqrt(HadronMassT2+Pz2)
602 // G4double zMin = (std::sqrt(HadronMa
603 G4double zMax = (std::sqrt(HadronMassT2+Pz2)
604
605 if (zMin >= zMax) return 0; // have to sta
606
607 G4double z = GetLightConeZ(zMin, zMax,
608 string->GetDecayParton()->Get
609 HadronPt.x(), HadronPt.y());
610
611 //... now compute hadron longitudinal moment
612 // longitudinal hadron momentum component Ha
613
614 HadronPt.setZ(0.5* string->GetDecayDirection
615 (z * string->LightConeDecay() -
616 HadronMassT2/(z * string->LightConeD
617 G4double HadronE = 0.5* (z * string->LightC
618 HadronMassT2/(z * string->LightConeD
619
620 G4LorentzVector * a4Momentum= new G4LorentzV
621
622 #ifdef debug_LUNDfragmentation
623 G4cout<<G4endl<<" string->GetDecayDire
624 G4cout<<"string->LightConeDecay() "<<s
625 G4cout<<"HadronPt,HadronE "<<HadronPt<
626 G4cout<<"String4Momentum "<<String4Mom
627 G4cout<<"Out of LUND SplitEandP "<<G4e
628 #endif
629
630 return a4Momentum;
631 }
632
633 //--------------------------------------------
634
635 G4double G4LundStringFragmentation::GetLightCo
636 G4int PD
637 G4Partic
638 G4double
639 {
640 G4double Mass = pHadron->GetPDGMass();
641 G4int HadronEncoding=std::abs(pHadron-
642
643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
644
645 G4double Alund, Blund;
646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(
647
648 if (!((std::abs(PDGEncodingOfDecayParton) >
649 { // ---------------- Quark fragmentation
650 Alund=1.;
651 Blund=0.7/GeV/GeV;
652
653 G4double BMt2 = Blund*Mt2;
654 if (Alund == 1.0) {
655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
656 else {
657 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sq
658 }
659
660 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
661 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
662 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blun
663
664 const G4int maxNumberOfLoops = 1000
665 G4int loopCounter = 0;
666 do
667 {
668 z = zmin + G4UniformRand()*(zmax-zmin);
669 //yf = (1-z)/z * G4Exp(-Blund*
670 yf = G4Pow::GetInstance()->powA(1.0-z,Alun
671 }
672 while ( (G4UniformRand()*maxYf > yf) && +
673 if ( loopCounter >= maxNumberOfLoop
674 z = 0.5*(zmin + zmax); // Just a
675 }
676 return z;
677 }
678
679 if (std::abs(PDGEncodingOfDecayParton) > 100
680 {
681 G4double an = 2.5;
682 an +=(sqr(Px)+sqr(Py))/sqr(GeV
683 z=zmin + (zmax-zmin)*G4Pow::Ge
684 if( PDGEncodingOfDecayParton >
685 }
686
687 return z;
688 }
689
690 //--------------------------------------------
691
692 G4bool G4LundStringFragmentation::SplitLast(G4
693 G4
694 G4
695 {
696 //... perform last cluster decay
697 SetMinimalStringMass( string);
698 if ( MinimalStringMass < 0.) return fa
699 #ifdef debug_LUNDfragmentation
700 G4cout<<G4endl<<"Split last-----------
701 G4cout<<"MinimalStringMass "<<MinimalS
702 G4cout<<"Left "<<string->GetLeftParto
703 G4cout<<"Right "<<string->GetRightPart
704 G4cout<<"String4mom "<<string->GetPstr
705 #endif
706
707 G4LorentzVector Str4Mom=string->Get4Mo
708 G4LorentzRotation toCms(-1*Str4Mom.boo
709 G4LorentzVector Pleft = toCms * string
710 toCms.rotateZ(-1*Pleft.phi());
711 toCms.rotateY(-1*Pleft.theta());
712
713 G4LorentzRotation toObserverFrame= toC
714
715 G4double StringMass=string->Mass();
716
717 G4ParticleDefinition * LeftHadron(0), * Righ
718
719 NumberOf_FS=0;
720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0.
721
722 G4int sampledState = 0;
723
724 #ifdef debug_LUNDfragmentation
725 G4cout<<"StrMass "<<StringMass<<" q's
726 <<string->GetLeftParton()->GetPa
727 <<string->GetRightParton()->GetP
728 #endif
729
730 string->SetLeftPartonStable(); // to query q
731
732 if (string->IsAFourQuarkString() )
733 {
734 G4int IDleft =std::abs(string->GetLe
735 G4int IDright=std::abs(string->GetRi
736
737 if ( (IDleft > 3000) || (IDright > 3
738 if ( ! Diquark_AntiDiquark_belowTh
739 {
740 return false;
741 }
742 } else {
743 // The string is qq-qqbar type. Diquarks a
744 if (StringMass-MinimalStringMass < 0
745 {
746 if (! Diquark_AntiDiquark_belowThreshold
747 {
748 return false;
749 }
750 } else
751 {
752 Diquark_AntiDiquark_aboveThreshold_lastS
753
754 if (NumberOf_FS == 0) return false;
755
756 sampledState = SampleS
757 if (string->GetLeftParton()->GetPDGEncod
758 {
759 LeftHadron =FS_LeftHadron[sampledState
760 RightHadron=FS_RightHadron[sampledStat
761 } else
762 {
763 LeftHadron =FS_RightHadron[sampledStat
764 RightHadron=FS_LeftHadron[sampledState
765 }
766 }
767 } // ID > 3300
768 } else {
769 if (string->DecayIsQuark() && string->Stab
770 { //... there are quarks on cluster
771 #ifdef debug_LUNDfragm
772 G4cout<<"Q Q string La
773 #endif
774
775 Quark_AntiQuark_lastSplitting(string, Le
776
777 if (NumberOf_FS == 0) return false;
778 sampledState = SampleState()
779 if (string->GetLeftParton()->GetPDGEncod
780 {
781 LeftHadron =FS_RightHadron[sampledStat
782 RightHadron=FS_LeftHadron[sampledState
783 } else
784 {
785 LeftHadron =FS_LeftHadron[sampledState
786 RightHadron=FS_RightHadron[sampledStat
787 }
788 } else
789 { //... there is a Diquark on one of
790 #ifdef debug_LUNDfragm
791 G4cout<<"DiQ Q string
792 #endif
793
794 Quark_Diquark_lastSplitting(string, Left
795
796 if (NumberOf_FS == 0) return false;
797 sampledState = SampleState()
798
799 if (string->GetLeftParton()->GetParticle
800 {
801 LeftHadron =FS_LeftHadron[sampledState
802 RightHadron=FS_RightHadron[sampledStat
803 } else
804 {
805 LeftHadron =FS_RightHadron[sampledStat
806 RightHadron=FS_LeftHadron[sampledState
807 }
808 }
809
810 }
811
812 #ifdef debug_LUNDfragmentation
813 G4cout<<"Sampled hadrons: "<<LeftHadro
814 #endif
815
816 G4LorentzVector P_left =string->GetPleft(),
817
818 G4LorentzVector LeftMom, RightMom;
819 G4ThreeVector Pos;
820
821 Sample4Momentum(&LeftMom, LeftHadron->GetPD
822 &RightMom, RightHadron->GetPDGMass(),
823 StringMass);
824
825 // Sample4Momentum ascribes LeftMom.pz
826 // It must be negative in case when th
827
828 if (!(string->DecayIsQuark() && string->Stab
829 { // Only for qq - q, q - qq, and qq - qqbar
830
831 if ( G4UniformRand() <= 0.5 )
832 {
833 if (P_left.z() <= 0.) {G4LorentzVector t
834 }
835 else
836 {
837 if (P_right.z() >= 0.) {G4LorentzVector
838 }
839 }
840
841 LeftMom *=toObserverFrame;
842 RightMom*=toObserverFrame;
843
844 LeftVector->push_back(new G4KineticTrack(Lef
845 RightVector->push_back(new G4KineticTrack(Ri
846
847 string->LorentzRotate(toObserverFrame);
848 return true;
849 }
850
851 //--------------------------------------------
852
853 G4bool G4LundStringFragmentation::
854 Diquark_AntiDiquark_belowThreshold_lastSplitti
855
856
857 {
858 G4double StringMass = string->Mass();
859
860 G4int cClusterInterrupt = 0;
861 G4bool isOK = false;
862 do
863 {
864 G4int LeftQuark1= string->GetLeftParton()-
865 G4int LeftQuark2=(string->GetLeftParton()-
866
867 G4int RightQuark1= string->GetRightParton(
868 G4int RightQuark2=(string->GetRightParton(
869
870 if (G4UniformRand()<0.5)
871 {
872 LeftHadron =hadronizer->Build(FindPartic
873 FindParticle(RightQuark1));
874 RightHadron= (LeftHadron == nullptr) ? n
875
876 FindParticle(RightQuark2));
877 } else
878 {
879 LeftHadron =hadronizer->Build(FindPartic
880 FindParticle(RightQuark2));
881 RightHadron=(LeftHadron == nullptr) ? nu
882 hadronizer
883 FindParticle(RightQuark1));
884 }
885
886 isOK = (LeftHadron != nullptr) && (RightHa
887
888 if(isOK) { isOK = (StringMass > LeftHadron
889 ++cClusterInterrupt;
890 //... repeat procedure, if mass of cluster
891 //... ClusterMassCut = 0.15*GeV model para
892 }
893 while (isOK == false && cClusterInterrupt <
894 /* Loop checking, 07.08.2015, A.Ribon */
895 return isOK;
896 }
897
898 //--------------------------------------------
899
900 G4bool G4LundStringFragmentation::
901 Diquark_AntiDiquark_aboveThreshold_lastSplitti
902
903
904 {
905 // StringMass-MinimalStringMass > 0. Creatio
906
907 G4double StringMass = string->Mass();
908 G4double StringMassSqr= sqr(StringMass);
909 G4ParticleDefinition * Di_Quark;
910 G4ParticleDefinition * Anti_Di_Quark;
911
912 if (string->GetLeftParton()->GetPDGEncoding(
913 {
914 Anti_Di_Quark =string->GetLeftParton();
915 Di_Quark=string->GetRightParton();
916 } else
917 {
918 Anti_Di_Quark =string->GetRightParton();
919 Di_Quark=string->GetLeftParton();
920 }
921
922 G4int IDAnti_di_quark =Anti_Di_Quark->Get
923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di
924 G4int IDdi_quark =Di_Quark->GetPDGEn
925 G4int AbsIDdi_quark =std::abs(IDdi_quar
926
927 G4int ADi_q1=AbsIDAnti_di_quark/1000;
928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000
929
930 G4int Di_q1=AbsIDdi_quark/1000;
931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
932
933 NumberOf_FS=0;
934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++)
935 {
936 G4int StateADiQ=0;
937 const G4int maxNumberOfLoops =
938 G4int loopCounter = 0;
939 do // while(Meson[AbsIDquark-1][ProdQ-1][
940 {
941 LeftHadron=G4ParticleTable::GetParticleT
942 -Baryon[ADi_q1-1][ADi_q2-1][Prod
943
944 if (LeftHadron == NULL) continue;
945 G4double LeftHadronMass=LeftHadron->GetP
946
947 G4int StateDiQ=0;
948 const G4int maxNumberO
949 G4int internalLoopCoun
950 do // while(Baryon[Di_q1-1][Di_q2-1][Pro
951 {
952 RightHadron=G4ParticleTable::GetPartic
953 +Baryon[Di_q1-1][Di_q2-1][Prod
954
955 if (RightHadron == NULL) continue;
956 G4double RightHadronMass=RightHadron->
957
958 if (StringMass > LeftHadronMass + Righ
959 {
960 if ( N
961 G4Ex
962 ed <
963 G4Ex
964
965 Numb
966 }
967
968 G4double FS_Psqr=lambda(StringMassSq
969 sqr(RightHadronMass));
970 //FS_Psqr=1.;
971 FS_Weight[NumberOf_FS]=std::sqrt(FS_
972 BaryonWeight[ADi_q1-1][AD
973 BaryonWeight[Di_q1-1][Di_
974 Prob_QQbar[ProdQ-1];
975
976 FS_LeftHadron[NumberOf_FS] = LeftHad
977 FS_RightHadron[NumberOf_FS]= RightHa
978 NumberOf_FS++;
979 } // End of if (StringMass > LeftHadro
980
981 StateDiQ++;
982
983 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ
984 ++internalLoo
985 if ( internalLoopCount
986 return false;
987 }
988
989 StateADiQ++;
990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ
991 ++loopCounter < maxNu
992 if ( loopCounter >= maxNumberO
993 return false;
994 }
995 } // End of for (G4int ProdQ=1; ProdQ < 4; P
996
997 return true;
998 }
999
1000 //-------------------------------------------
1001
1002 G4bool G4LundStringFragmentation::Quark_Diqua
1003
1004
1005 {
1006 G4double StringMass = string->Mass();
1007 G4double StringMassSqr= sqr(StringMass);
1008
1009 G4ParticleDefinition * Di_Quark;
1010 G4ParticleDefinition * Quark;
1011
1012 if (string->GetLeftParton()->GetParticleSub
1013 {
1014 Quark =string->GetLeftParton();
1015 Di_Quark=string->GetRightParton();
1016 } else
1017 {
1018 Quark =string->GetRightParton();
1019 Di_Quark=string->GetLeftParton();
1020 }
1021
1022 G4int IDquark =Quark->GetPDGEncoding
1023 G4int AbsIDquark =std::abs(IDquark);
1024 G4int IDdi_quark =Di_Quark->GetPDGEncodin
1025 G4int AbsIDdi_quark=std::abs(IDdi_quark);
1026 G4int Di_q1=AbsIDdi_quark/1000;
1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1028 G4int SignDiQ= 1;
1029 if (IDdi_quark < 0) SignDiQ=-1;
1030
1031 NumberOf_FS=0;
1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) //
1033 { //
1034 G4int SignQ;
1035 if (IDquark > 0)
1036 {
1037 SignQ=-1;
1038 if (IDquark == 2)
1039 if ((IDquark == 1) && (ProdQ == 3)) Sig
1040 if ((IDquark == 3) && (ProdQ == 1)) Sig
1041 if (IDquark == 4)
1042 if (IDquark == 5) Sig
1043 } else
1044 {
1045 SignQ= 1;
1046 if (IDquark == -2) Sig
1047 if ((IDquark ==-1) && (ProdQ == 3)) Sig
1048 if ((IDquark ==-3) && (ProdQ == 1)) Sig
1049 if (IDquark == -4) Sig
1050 if (IDquark == -5) Sig
1051 }
1052
1053 if (AbsIDquark == ProdQ) SignQ
1054
1055 G4int StateQ=0;
1056 const G4int maxNumberOfLoops
1057 G4int loopCounter = 0;
1058 do // while(Meson[AbsIDquark-1][ProdQ-1]
1059 {
1060 LeftHadron=G4ParticleTable::GetParticle
1061 Meson[AbsIDquark-1][ProdQ-1][St
1062 if (LeftHadron == NULL) continue;
1063 G4double LeftHadronMass=LeftHadron->Get
1064
1065 G4int StateDiQ=0;
1066 const G4int maxNumber
1067 G4int internalLoopCou
1068 do // while(Baryon[Di_q1-1][Di_q2-1][Pr
1069 {
1070 RightHadron=G4ParticleTable::GetParti
1071 Baryon[Di_q1-1][Di_q2-1][Prod
1072 if (RightHadron == NULL) continue;
1073 G4double RightHadronMass=RightHadron-
1074
1075 if (StringMass > LeftHadronMass + Rig
1076 {
1077 if (
1078 G4E
1079 ed
1080 G4E
1081
1082 Num
1083 }
1084
1085 G4double FS_Psqr=lambda(StringMassS
1086 sqr(RightHadronMass));
1087 FS_Weight[NumberOf_FS]=std::sqrt(FS
1088 MesonWeight[AbsIDquark-1
1089 BaryonWeight[Di_q1-1][Di
1090 Prob_QQbar[ProdQ-1];
1091
1092 FS_LeftHadron[NumberOf_FS] = LeftHa
1093 FS_RightHadron[NumberOf_FS]= RightH
1094
1095 NumberOf_FS++;
1096 } // End of if (StringMass > LeftHadr
1097
1098 StateDiQ++;
1099
1100 } while( (Baryon[Di_q1-1][Di_q2-1][Prod
1101 ++internalLo
1102 if ( internalLoopCoun
1103 return false;
1104 }
1105
1106 StateQ++;
1107 } while( (Meson[AbsIDquark-1][ProdQ-1][St
1108 ++loopCounter < maxN
1109
1110 if ( loopCounter >= maxNumb
1111 return false;
1112 }
1113 }
1114
1115 return true;
1116 }
1117
1118 //-------------------------------------------
1119
1120 G4bool G4LundStringFragmentation::Quark_AntiQ
1121
1122
1123 {
1124 G4double StringMass = string->Mass();
1125 G4double StringMassSqr= sqr(StringMass);
1126
1127 G4ParticleDefinition * Quark;
1128 G4ParticleDefinition * Anti_Quark;
1129
1130 if (string->GetLeftParton()->GetPDGEncoding
1131 {
1132 Quark =string->GetLeftParton();
1133 Anti_Quark=string->GetRightParton();
1134 } else
1135 {
1136 Quark =string->GetRightParton();
1137 Anti_Quark=string->GetLeftParton();
1138 }
1139
1140 G4int IDquark =Quark->GetPDGEncodin
1141 G4int AbsIDquark =std::abs(IDquark);
1142 G4int QuarkCharge =Qcharge[IDquar
1143
1144 G4int IDanti_quark =Anti_Quark->GetPDGEn
1145 G4int AbsIDanti_quark =std::abs(IDanti_quar
1146 G4int AntiQuarkCharge =-Qcharge[AbsID
1147
1148 G4int LeftHadronCharge(0), RightHadro
1149
1150 //G4cout<<"Q Qbar "<<IDquark<<" "<<ID
1151
1152 NumberOf_FS=0;
1153 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) //
1154 { //
1155 //G4cout<<"NumberOf_FS ProdQ
1156 LeftHadronCharge = QuarkCharge - Qcharge[
1157 G4int SignQ = LeftHadronCharge/3; if (Sig
1158
1159 if ((IDquark == 1) && (ProdQ == 3)) SignQ
1160 if ((IDquark == 3) && (ProdQ == 1)) SignQ
1161 if ((IDquark == 4) && (ProdQ
1162 if ((IDquark == 5) && (ProdQ
1163 if ((IDquark == 5) && (ProdQ
1164
1165 RightHadronCharge = AntiQuark
1166 G4int SignAQ = RightHadronCharge/3; if (S
1167
1168 if ((IDanti_quark ==-1) && (ProdQ == 3))
1169 if ((IDanti_quark ==-3) && (ProdQ == 1))
1170 if ((IDanti_quark ==-4) && (ProdQ == 2))
1171 if ((IDanti_quark ==-5) && (ProdQ == 1))
1172 if ((IDanti_quark ==-5) && (ProdQ == 3))
1173
1174 //G4cout<<"ProQ signs "<<Prod
1175
1176 G4int StateQ=0;
1177 const G4int maxNumberOfLoops
1178 G4int loopCounter = 0;
1179 do
1180 {
1181 //G4cout<<"[AbsIDquar
1182 //<<ProdQ-1<<" "<<Sta
1183 LeftHadron=G4ParticleTable::GetParticle
1184 Meson[AbsIDquark-1][ProdQ-
1185 //G4cout<<"LeftHadron
1186 if (LeftHadron == NULL) { StateQ++; con
1187 //G4cout<<"LeftHadron
1188 G4double LeftHadronMass=LeftHadron->Get
1189
1190 G4int StateAQ=0;
1191 const G4int maxNumber
1192 G4int internalLoopCou
1193 do
1194 {
1195 //G4cout<<"
1196 // <<Pro
1197 RightHadron=G4ParticleTable::GetParti
1198 Meson[AbsIDanti_quark-1][Prod
1199 //G4cout<<"Ri
1200 if(RightHadron == NULL) { StateAQ++;
1201 //G4cout<<"Ri
1202 G4double RightHadronMass=RightHadron-
1203
1204 if (StringMass > LeftHadronMass + Rig
1205 {
1206 if (
1207 G4E
1208 ed
1209 G4E
1210
1211 Num
1212 }
1213
1214 G4dou
1215 sqr(RightHadronMass));
1216 //FS_Psqr=1.;
1217 FS_Weight[NumberOf_FS]=std::sqrt(FS
1218 MesonWeight[AbsIDquark-1
1219 MesonWeight[AbsIDanti_qu
1220 Prob_QQbar[ProdQ-1];
1221 if (string->GetLeftParton()->GetPDG
1222 {
1223 FS_LeftHadron[NumberOf_FS] = Righ
1224 FS_RightHadron[NumberOf_FS]= Left
1225 } else
1226 {
1227 FS_LeftHadron[NumberOf_FS] = Left
1228 FS_RightHadron[NumberOf_FS]= Righ
1229 }
1230
1231 NumberOf_FS++;
1232
1233 }
1234
1235 StateAQ++;
1236 //G4cout<<"
1237 // <<Mes
1238 } while ( (Meson[AbsIDanti_quark-1][Pro
1239 ++internalL
1240 if ( internalLoopCo
1241 return false;
1242 }
1243
1244 StateQ++;
1245 //G4cout<<"StateQ Mes
1246 // <<Meson[AbsID
1247
1248 } while ( (Meson[AbsIDquark-1][ProdQ-1][S
1249 ++loopCounter < maxN
1250 if ( loopCounter >= maxNumb
1251 return false;
1252 }
1253 } // End of for (G4int ProdQ=1; ProdQ < 4;
1254
1255 return true;
1256 }
1257
1258 //-------------------------------------------
1259
1260 G4int G4LundStringFragmentation::SampleState(
1261 {
1262 if ( NumberOf_FS > 349 ) {
1263 G4ExceptionDescription ed;
1264 ed << " NumberOf_FS exceeds its lim
1265 G4Exception( "G4LundStringFragmenta
1266 NumberOf_FS = 349;
1267 }
1268
1269 G4double SumWeights=0.;
1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWei
1271
1272 G4double ksi=G4UniformRand();
1273 G4double Sum=0.;
1274 G4int indexPosition = 0;
1275
1276 for (G4int i=0; i<NumberOf_FS; i++)
1277 {
1278 Sum+=(FS_Weight[i]/SumWeights);
1279 indexPosition=i;
1280 if (Sum >= ksi) break;
1281 }
1282 return indexPosition;
1283 }
1284
1285 //-------------------------------------------
1286
1287 void G4LundStringFragmentation::Sample4Moment
1288
1289
1290 {
1291 // ------ Sampling of momenta of 2 last pro
1292 G4ThreeVector Pt;
1293 G4double MassMt, AntiMassMt;
1294 G4double AvailablePz, AvailablePz2;
1295 #ifdef debug_LUNDfragmentation
1296 G4cout<<"Sampling of momenta of 2 las
1297 G4cout<<"Init Mass "<<InitialMass<<"
1298 #endif
1299
1300 G4double r_val = sqr(InitialMass*InitialMas
1301 sqr(2.*Mass*AntiMass);
1302 G4double Pabs = (r_val > 0.)? std::sqrt(r_v
1303
1304 const G4int maxNumberOfLoops = 1000;
1305 G4double SigmaQTw = SigmaQT;
1306 if ( Mass > 930. || AntiMass > 930. ) {
1307 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMa
1308 }
1309 if ( Mass < 930. && AntiMass < 930. )
1310 if ( ( Mass < 930. && AntiMass > 930.
1311 ( Mass > 930. && AntiMass < 930. ) ) {
1312 //SigmaQT = -1.; // isotropical de
1313 }
1314 if ( Mass > 930. && AntiMass > 930. )
1315 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+
1316 }
1317
1318 G4int loopCounter = 0;
1319 do
1320 {
1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou
1322 MassMt = std::sqrt( Mass * Mass
1323 AntiMassMt= std::sqrt(AntiMass * AntiMass
1324 }
1325 while ( (InitialMass < MassMt + AntiMassMt)
1326
1327 SigmaQT = SigmaQTw;
1328
1329 AvailablePz2= sqr(InitialMass*InitialMass -
1330 4.*sqr(MassMt*AntiMassMt);
1331
1332 AvailablePz2 /=(4.*InitialMass*InitialMass)
1333 AvailablePz = std::sqrt(AvailablePz2);
1334
1335 G4double Px=Pt.getX();
1336 G4double Py=Pt.getY();
1337
1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(
1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz
1340
1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A
1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av
1343
1344 #ifdef debug_LUNDfragmentation
1345 G4cout<<"Fmass Mom "<<Mom->getX()<<"
1346 G4cout<<"Smass Mom "<<AntiMom->getX()
1347 <<" "<<AntiMom->getT()<<G4endl;
1348 #endif
1349 }
1350
1351 //-------------------------------------------
1352
1353 G4double G4LundStringFragmentation::lambda(G4
1354 {
1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4
1356 return lam;
1357 }
1358
1359 // ------------------------------------------
1360 G4LundStringFragmentation::~G4LundStringFragm
1361 {}
1362
1363