Geant4-11
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//
28// -----------------------------------------------------------------------------
29// GEANT 4 class implementation file
30//
31// History: first implementation, Maxim Komogorov, 10-Jul-1998
32// -----------------------------------------------------------------------------
35#include "G4SystemOfUnits.hh"
36#include "Randomize.hh"
38#include "G4DiQuarks.hh"
39#include "G4Quarks.hh"
41#include "G4Exp.hh"
42#include "G4Pow.hh"
43
44//#define debug_LUNDfragmentation
45
46// Class G4LundStringFragmentation
47//*************************************************************************************
48
50 : G4VLongitudinalStringDecay("LundStringFragmentation")
51{
52 SetMassCut(210.*MeV); // Mpi + Delta
53 // For ProduceOneHadron it is required
54 // that no one pi-meson can be produced.
55 SigmaQT = 0.435 * GeV;
56 Tmt = 190.0 * MeV;
59 SetStrangenessSuppression((1.0 - 0.12)/2.0);
61
62 // Check if charmed and bottom hadrons are enabled: if this is the case, then
63 // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum,
64 // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom
65 // hadrons can't/can be created during the string fragmentation of ordinary
66 // (i.e. not heavy) projectile hadron nuclear reactions.
67 if ( G4HadronicParameters::Instance()->EnableBCParticles() ) {
68 SetProbCCbar(0.005); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
69 SetProbBBbar(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
70 } else {
71 SetProbCCbar(0.0);
72 SetProbBBbar(0.0);
73 }
74
75 SetMinMasses(); // For treating of small string decays
76}
77
78//--------------------------------------------------------------------------------------
79
81{
82 // Can no longer modify Parameters for Fragmentation.
83
84 PastInitPhase=true;
85
86 G4FragmentingString aString(theString);
87 SetMinimalStringMass(&aString);
88
89 #ifdef debug_LUNDfragmentation
90 G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl;
91 G4cout<<G4endl<<"LUND StringFragm: String Mass "
92 <<theString.Get4Momentum().mag()<<G4endl
93 <<"4Mom "<<theString.Get4Momentum()<<G4endl
94 <<"------------------------------------"<<G4endl;
95 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
96 <<theString.GetRightParton()->GetPDGcode()<<" "
97 <<theString.GetDirection()<< G4endl;
98 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
99 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl<<G4endl;
100 G4cout<<"Check for Fragmentation "<<G4endl;
101 #endif
102
103 G4KineticTrackVector * LeftVector(0);
104
105 if (!aString.IsAFourQuarkString() && !IsItFragmentable(&aString))
106 {
107 #ifdef debug_LUNDfragmentation
108 G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
109 #endif
110 // SetMassCut(210.*MeV); // For ProduceOneHadron it is required
111 // that no one pi-meson can be produced.
112
113 G4double Mcut = GetMassCut();
114 SetMassCut(10000.*MeV);
115 LeftVector=ProduceOneHadron(&theString);
116 SetMassCut(Mcut);
117
118 if ( LeftVector )
119 {
120 if ( LeftVector->size() > 0)
121 {
122 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
123 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
124 }
125 if (LeftVector->size() > 1)
126 {
127 // 2 hadrons created from qq-qqbar are stored
128 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
129 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
130 }
131 }
132 return LeftVector;
133 }
134
135 #ifdef debug_LUNDfragmentation
136 G4cout<<"The string will be fragmented. "<<G4endl;
137 #endif
138
139 // The string can fragment. At least two particles can be produced.
140 LeftVector =new G4KineticTrackVector;
142
143 G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector);
144
145 if ( ! success )
146 {
147 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
148 LeftVector->clear();
149 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
150 delete RightVector;
151 return LeftVector;
152 }
153
154 // Join Left- and RightVector into LeftVector in correct order.
155 while (!RightVector->empty())
156 {
157 LeftVector->push_back(RightVector->back());
158 RightVector->erase(RightVector->end()-1);
159 }
160 delete RightVector;
161
162 return LeftVector;
163}
164
165//----------------------------------------------------------------------------------
166
168{
169 SetMinimalStringMass(string);
170 //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
171
172 return std::abs(MinimalStringMass) < string->Get4Momentum().mag();
173
174 //MinimalStringMass is negative and large for a string with unknown particles in a final 2-particle decay.
175}
176
177//----------------------------------------------------------------------------------------
178
180 G4KineticTrackVector * & LeftVector,
181 G4KineticTrackVector * & RightVector )
182{
183 #ifdef debug_LUNDfragmentation
184 G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" "
185 <<theString.GetLeftParton()->Get4Momentum()<<G4endl
186 <<" "<<theString.GetRightParton()->GetPDGcode()<<" "
187 <<theString.GetRightParton()->Get4Momentum()<<G4endl
188 <<"Direction "<<theString.GetDirection()<< G4endl;
189 #endif
190
191 G4LorentzRotation toCmsI, toObserverFrameI;
192
193 G4bool final_success=false;
194 G4bool inner_success=true;
195
196 G4int attempt=0;
197
198 while ( ! final_success && attempt++ < StringLoopInterrupt )
199 { // If the string fragmentation does not be happend,
200 // repeat the fragmentation.
201
202 G4FragmentingString *currentString = new G4FragmentingString(theString);
203 toCmsI = currentString->TransformToAlignedCms();
204 toObserverFrameI = toCmsI.inverse();
205
206 G4LorentzRotation toCms, toObserverFrame;
207
208 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
209
210 // Cleaning up the previously produced hadrons
211 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
212 LeftVector->clear();
213 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
214 RightVector->clear();
215
216 // Main fragmentation loop until the string will not be able to fragment
217 inner_success=true; // set false on failure.
218 const G4int maxNumberOfLoops = 1000;
219 G4int loopCounter = -1;
220
221 while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
222 { // Split current string into hadron + new string
223 #ifdef debug_LUNDfragmentation
224 G4cout<<"The string will fragment. "<<G4endl;;
225 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
226 #endif
227 G4FragmentingString *newString=0; // used as output from SplitUp.
228
229 toCms=currentString->TransformToAlignedCms();
230 toObserverFrame= toCms.inverse();
231
232 #ifdef debug_LUNDfragmentation
233 //G4cout<<"CMS Left mom "<<currentString->GetPleft()<<G4endl;
234 //G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl;
235 //G4cout<<"CMS String M "<<currentString->GetPstring()<<G4endl;
236 #endif
237
238 G4KineticTrack * Hadron=Splitup(currentString,newString);
239
240 if ( Hadron != 0 ) // Store the hadron
241 {
242 #ifdef debug_LUNDfragmentation
243 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
244 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
245 #endif
246
247 Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum());
248
249 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
250 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
251
252 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
253 G4LorentzVector Momentum = toObserverFrame*Coordinate;
254 Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
255 G4ThreeVector aPosition(Momentum.vect());
256 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
257
258 // Open to protect hadron production at fragmentation
259 if ( currentString->GetDecayDirection() > 0 )
260 {
261 LeftVector->push_back(Hadron);
262 } else
263 {
264 RightVector->push_back(Hadron);
265 }
266 delete currentString;
267 currentString=newString;
268 } else {
269 if ( newString ) delete newString;
270 }
271
272 currentString->LorentzRotate(toObserverFrame);
273 };
274
275 if ( loopCounter >= maxNumberOfLoops ) {
276 inner_success=false;
277 }
278
279 // Split remaining string into 2 final hadrons.
280 #ifdef debug_LUNDfragmentation
281 if (inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
282 #endif
283
284 if ( inner_success && SplitLast(currentString, LeftVector, RightVector) ) // Close to protect Last Str. Decay
285 {
286 final_success = true;
287 }
288
289 delete currentString;
290 } // End of the loop where we try to fragment the string.
291
292 G4int sign = +1;
293 if ( theString.GetDirection() < 0 ) sign = -1;
294 for ( unsigned int hadronI = 0; hadronI < LeftVector->size(); ++hadronI ) {
295 G4LorentzVector Tmp = LeftVector->operator[](hadronI)->Get4Momentum();
296 Tmp.setZ(sign*Tmp.getZ());
297 Tmp *= toObserverFrameI;
298 LeftVector->operator[](hadronI)->Set4Momentum(Tmp);
299 }
300 for ( unsigned int hadronI = 0; hadronI < RightVector->size(); ++hadronI ) {
301 G4LorentzVector Tmp = RightVector->operator[](hadronI)->Get4Momentum();
302 Tmp.setZ(sign*Tmp.getZ());
303 Tmp *= toObserverFrameI;
304 RightVector->operator[](hadronI)->Set4Momentum(Tmp);
305 }
306
307 return final_success;
308}
309
310//----------------------------------------------------------------------------------------
311
313{
314 SetMinimalStringMass(string);
315
316 if ( MinimalStringMass < 0.) return true;
317
318 if (string->IsAFourQuarkString())
319 {
320 return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
321 } else {
322
323 if (MinimalStringMass < 0.0 ) return false; // For a string with di-quark having c or b quarks and s, c, b quarks
324
325 G4bool Result = G4UniformRand() <
326 G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
327 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND
328
329 #ifdef debug_LUNDfragmentation
330 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass
331 <<" "<<string->Mass()<<G4endl;
332 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
333 #endif
334 return Result;
335 }
336}
337
338//-----------------------------------------------------------------------------
339
341 G4FragmentingString *&newString)
342{
343 #ifdef debug_LUNDfragmentation
344 G4cout<<G4endl;
345 G4cout<<"Start SplitUP ========================="<<G4endl;
346 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
347 <<string->GetRightParton()->GetPDGEncoding()<<" "
348 <<"Direction " <<string->GetDecayDirection()<<G4endl;
349 #endif
350
351 //... random choice of string end to use for creating the hadron (decay)
352 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
353 if (SideOfDecay < 0)
354 {
355 string->SetLeftPartonStable();
356 } else
357 {
358 string->SetRightPartonStable();
359 }
360
361 G4ParticleDefinition *newStringEnd;
362 G4ParticleDefinition * HadronDefinition;
363
364 G4double StringMass=string->Mass();
365
366 G4double ProbDqADq = GetDiquarkSuppress();
367 G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
368
369 #ifdef debug_LUNDfragmentation
370 G4cout<<"StrMass DiquarkSuppression "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
371 #endif
372
373 G4int NumberOfpossibleBaryons = 2;
374
375 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
376 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
377
378 G4double ActualProb = ProbDqADq ;
379 ActualProb *= (1.0-sqr(NumberOfpossibleBaryons*1400.0/StringMass));
380
381 SetDiquarkSuppression(ActualProb);
382
383 G4double Mth = 1250.0; // 2 Mk + Mpi
384 if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;} // Mlambda/Msigma + Mk + Mpi
385 else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;} // 2 Mlambda/Msigma + Mk + Mpi
386 else {}
387
388 ActualProb = ProbSaS * (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass , 4.0 ));
389 if ( ActualProb < 0.0 ) ActualProb = 0.0;
390 SetStrangenessSuppression((1.0-ActualProb)/2.0);
391
392 #ifdef debug_LUNDfragmentation
393 G4cout<<"StrMass DiquarkSuppression corrected "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
394 #endif
395
396 if (string->DecayIsQuark())
397 {
398 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
399 } else {
400 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
401 }
402
403 SetDiquarkSuppression(ProbDqADq);
404 SetStrangenessSuppression((1.0-ProbSaS)/2.0);
405
406 if ( HadronDefinition == NULL ) { G4KineticTrack * Hadron =0; return Hadron; }
407
408 #ifdef debug_LUNDfragmentation
409 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
410 <<" produces hadron "<<HadronDefinition->GetParticleName()
411 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
412 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
413 #endif
414 // create new String from old, ie. keep Left and Right order, but replace decay
415
416 if ( newString ) delete newString;
417
418 newString=new G4FragmentingString(*string,newStringEnd); // To store possible quark containt of new string
419
420 #ifdef debug_LUNDfragmentation
421 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
422 #endif
423 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
424
425 delete newString; newString=0;
426
427 G4KineticTrack * Hadron =0;
428 if ( HadronMomentum != 0 ) {
429
430 #ifdef debug_LUNDfragmentation
431 G4cout<<"The attempt was successful"<<G4endl;
432 #endif
434 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
435
436 if ( newString ) delete newString;
437
438 newString=new G4FragmentingString(*string,newStringEnd,
439 HadronMomentum);
440 delete HadronMomentum;
441 }
442 else
443 {
444 #ifdef debug_LUNDfragmentation
445 G4cout<<"The attempt was not successful !!!"<<G4endl;
446 #endif
447 }
448
449 #ifdef debug_LUNDfragmentation
450 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
451 #endif
452
453 return Hadron;
454}
455
456//-----------------------------------------------------------------------------
457
459 G4ParticleDefinition *&created)
460{
462 G4double ProbQQbar = (1.0 - 2.0*StrSup);
463
464 //... can Diquark break or not?
466
467 //... Diquark break
468 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
469 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
470 if (G4UniformRand() < 0.5)
471 {
472 G4int Swap = stableQuarkEncoding;
473 stableQuarkEncoding = decayQuarkEncoding;
474 decayQuarkEncoding = Swap;
475 }
476
477 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark
478
479 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
480
481 //... Build new Diquark
482 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
483 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
484 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
485 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
486 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
487 created = FindParticle(NewDecayEncoding);
488 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
489 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
490 StrangeSuppress=StrSup;
491
492 return had;
493
494 } else {
495 //... Diquark does not break
496
497 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark
498
499 StrangeSuppress=(1.0 - ProbQQbar * 0.9)/2.0;
500 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
501
502 created = QuarkPair.second;
503
504 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
505 StrangeSuppress=StrSup;
506
507 return had;
508 }
509}
510
511//-----------------------------------------------------------------------------
512
514 G4FragmentingString * string,
515 G4FragmentingString * newString)
516{
517 G4LorentzVector String4Momentum=string->Get4Momentum();
518 G4double StringMT2=string->MassT2();
519 G4double StringMT =std::sqrt(StringMT2);
520
521 G4double HadronMass = pHadron->GetPDGMass();
522 SetMinimalStringMass(newString);
523
524 if ( MinimalStringMass < 0.0 ) return nullptr;
525
526 #ifdef debug_LUNDfragmentation
527 G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
528 G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()
529 <<" "<<std::sqrt(StringMT2)<<G4endl;
530 G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
531 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
532 <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
533 #endif
534
535 if ((HadronMass + MinimalStringMass > string->Mass()) || MinimalStringMass < 0.)
536 {
537 #ifdef debug_LUNDfragmentation
538 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
539 #endif
540 return 0;
541 } // have to start all over!
542
543 String4Momentum.setPz(0.);
544 G4ThreeVector StringPt=String4Momentum.vect();
545 StringPt.setZ(0.);
546
547 // calculate and assign hadron transverse momentum component HadronPx and HadronPy
548 G4ThreeVector HadronPt , RemSysPt;
549 G4double HadronMassT2, ResidualMassT2;
550 G4double HadronMt, Pt, Pt2, phi;
551
552 G4double TmtCur = Tmt;
553 if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
554 (pHadron->GetBaryonNumber() != 0) ) {
555 TmtCur = Tmt*0.37; // q->B
556 } else if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
557 (pHadron->GetBaryonNumber() == 0) ) {
558 //TmtCur = Tmt; // q->M
559 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") &&
560 (pHadron->GetBaryonNumber() == 0) ) {
561 //TmtCur = Tmt*0.89; // qq -> M
562 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") &&
563 (pHadron->GetBaryonNumber() != 0) ) {
564 TmtCur = Tmt*1.35; // qq -> B
565 }
566
567 //... sample Pt of the hadron
568 G4int attempt=0;
569 do
570 {
571 attempt++; if (attempt > StringLoopInterrupt) {return 0;}
572
573 HadronMt = HadronMass - TmtCur*G4Log(G4UniformRand());
574 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
575 phi = 2.*pi*G4UniformRand();
576 HadronPt = G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0. );
577 RemSysPt = StringPt - HadronPt;
578 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
579 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
580
581 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
582
583 //... sample z to define hadron longitudinal momentum and energy
584 //... but first check the available phase space
585
586 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
587 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
588
589 if (Pz2 < 0 ) {return 0;} // have to start all over!
590
591 //... then compute allowed z region z_min <= z <= z_max
592
593 G4double Pz = std::sqrt(Pz2);
594 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
595 // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes
596 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
597
598 if (zMin >= zMax) return 0; // have to start all over!
599
600 G4double z = GetLightConeZ(zMin, zMax,
601 string->GetDecayParton()->GetPDGEncoding(), pHadron,
602 HadronPt.x(), HadronPt.y());
603
604 //... now compute hadron longitudinal momentum and energy
605 // longitudinal hadron momentum component HadronPz
606
607 HadronPt.setZ(0.5* string->GetDecayDirection() *
608 (z * string->LightConeDecay() -
609 HadronMassT2/(z * string->LightConeDecay())));
610 G4double HadronE = 0.5* (z * string->LightConeDecay() +
611 HadronMassT2/(z * string->LightConeDecay()));
612
613 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
614
615 #ifdef debug_LUNDfragmentation
616 G4cout<<G4endl<<" string->GetDecayDirection() "<<string->GetDecayDirection()<<G4endl<<G4endl;
617 G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
618 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
619 G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
620 G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl;
621 #endif
622
623 return a4Momentum;
624}
625
626//-----------------------------------------------------------------------------------------
627
629 G4int PDGEncodingOfDecayParton,
630 G4ParticleDefinition* pHadron,
631 G4double Px, G4double Py)
632{
633 G4double Mass = pHadron->GetPDGMass();
634 G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
635
636 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
637
638 G4double Alund, Blund;
639 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
640
641 if (!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
642 { // ---------------- Quark fragmentation and qq-> meson ----------------------
643 Alund=1.;
644 Blund=0.7/GeV/GeV;
645
646 G4double BMt2 = Blund*Mt2;
647 if (Alund == 1.0) {
648 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
649 else {
650 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
651 }
652
653 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
654 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
655 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf);
656
657 const G4int maxNumberOfLoops = 1000;
658 G4int loopCounter = 0;
659 do
660 {
661 z = zmin + G4UniformRand()*(zmax-zmin);
662 //yf = (1-z)/z * G4Exp(-Blund*Mt2/z);
663 yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z);
664 }
665 while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
666 if ( loopCounter >= maxNumberOfLoops ) {
667 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
668 }
669 return z;
670 }
671
672 if (std::abs(PDGEncodingOfDecayParton) > 1000)
673 {
674 G4double an = 2.5;
675 an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5;
676 z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an);
677 }
678
679 return z;
680}
681
682//----------------------------------------------------------------------------------------------------------
683
685 G4KineticTrackVector * LeftVector,
686 G4KineticTrackVector * RightVector)
687{
688 //... perform last cluster decay
689 SetMinimalStringMass( string);
690 if ( MinimalStringMass < 0.) return false;
691 #ifdef debug_LUNDfragmentation
692 G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl;
693 G4cout<<"MinimalStringMass "<<MinimalStringMass<<G4endl;
694 G4cout<<"Left "<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetPleft()<<G4endl;
695 G4cout<<"Right "<<string->GetRightParton()->GetPDGEncoding()<<" "<<string->GetPright()<<G4endl;
696 G4cout<<"String4mom "<<string->GetPstring()<<" "<<string->GetPstring().mag()<<G4endl;
697 #endif
698
699 G4LorentzVector Str4Mom=string->Get4Momentum();
700 G4LorentzRotation toCms(-1*Str4Mom.boostVector());
701 G4LorentzVector Pleft = toCms * string->GetPleft();
702 toCms.rotateZ(-1*Pleft.phi());
703 toCms.rotateY(-1*Pleft.theta());
704
705 G4LorentzRotation toObserverFrame= toCms.inverse();
706
707 G4double StringMass=string->Mass();
708
709 G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
710
711 NumberOf_FS=0;
712 for (G4int i=0; i<350; i++) {FS_Weight[i]=0.;}
713
714 G4int sampledState = 0;
715
716 #ifdef debug_LUNDfragmentation
717 G4cout<<"StrMass "<<StringMass<<" q's "
718 <<string->GetLeftParton()->GetParticleName()<<" "
719 <<string->GetRightParton()->GetParticleName()<<G4endl;
720 #endif
721
722 string->SetLeftPartonStable(); // to query quark contents..
723
724 if (string->IsAFourQuarkString() )
725 {
726 // The string is qq-qqbar type. Diquarks are on the string ends
727 if (StringMass-MinimalStringMass < 0.)
728 {
729 if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
730 {
731 return false;
732 }
733 } else
734 {
735 Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
736
737 if (NumberOf_FS == 0) return false;
738
739 sampledState = SampleState();
740 if (string->GetLeftParton()->GetPDGEncoding() < 0)
741 {
742 LeftHadron =FS_LeftHadron[sampledState];
743 RightHadron=FS_RightHadron[sampledState];
744 } else
745 {
746 LeftHadron =FS_RightHadron[sampledState];
747 RightHadron=FS_LeftHadron[sampledState];
748 }
749 }
750 } else
751 {
752 if (string->DecayIsQuark() && string->StableIsQuark() )
753 { //... there are quarks on cluster ends
754 #ifdef debug_LUNDfragmentation
755 G4cout<<"Q Q string LastSplit"<<G4endl;
756 #endif
757
758 Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
759
760 if (NumberOf_FS == 0) return false;
761 sampledState = SampleState();
762 if (string->GetLeftParton()->GetPDGEncoding() < 0)
763 {
764 LeftHadron =FS_RightHadron[sampledState];
765 RightHadron=FS_LeftHadron[sampledState];
766 } else
767 {
768 LeftHadron =FS_LeftHadron[sampledState];
769 RightHadron=FS_RightHadron[sampledState];
770 }
771 } else
772 { //... there is a Diquark on one of the cluster ends
773 #ifdef debug_LUNDfragmentation
774 G4cout<<"DiQ Q string Last Split"<<G4endl;
775 #endif
776
777 Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
778
779 if (NumberOf_FS == 0) return false;
780 sampledState = SampleState();
781
782 if (string->GetLeftParton()->GetParticleSubType() == "quark")
783 {
784 LeftHadron =FS_LeftHadron[sampledState];
785 RightHadron=FS_RightHadron[sampledState];
786 } else
787 {
788 LeftHadron =FS_RightHadron[sampledState];
789 RightHadron=FS_LeftHadron[sampledState];
790 }
791 }
792
793 }
794
795 #ifdef debug_LUNDfragmentation
796 G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
797 #endif
798
799 G4LorentzVector P_left =string->GetPleft(), P_right = string->GetPright();
800
801 G4LorentzVector LeftMom, RightMom;
803
804 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
805 &RightMom, RightHadron->GetPDGMass(),
806 StringMass);
807
808 // Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases.
809 // It must be negative in case when the system is moving against Z axis.
810
811 if (!(string->DecayIsQuark() && string->StableIsQuark() ))
812 { // Only for qq - q, q - qq, and qq - qqbar -------------------
813
814 if ( G4UniformRand() <= 0.5 )
815 {
816 if (P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
817 }
818 else
819 {
820 if (P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
821 }
822 }
823
824 LeftMom *=toObserverFrame;
825 RightMom*=toObserverFrame;
826
827 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
828 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
829
830 string->LorentzRotate(toObserverFrame);
831 return true;
832}
833
834//----------------------------------------------------------------------------------------
835
838 G4ParticleDefinition * & LeftHadron,
839 G4ParticleDefinition * & RightHadron)
840{
841 G4double StringMass = string->Mass();
842
843 G4int cClusterInterrupt = 0;
844 G4bool isOK = false;
845 do
846 {
847 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
848 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
849
850 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
851 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
852
853 if (G4UniformRand()<0.5)
854 {
855 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
856 FindParticle(RightQuark1));
857 RightHadron= (LeftHadron == nullptr) ? nullptr :
858 hadronizer->Build(FindParticle( LeftQuark2),
859 FindParticle(RightQuark2));
860 } else
861 {
862 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
863 FindParticle(RightQuark2));
864 RightHadron=(LeftHadron == nullptr) ? nullptr :
865 hadronizer->Build(FindParticle( LeftQuark2),
866 FindParticle(RightQuark1));
867 }
868
869 isOK = (LeftHadron != nullptr) && (RightHadron != nullptr);
870 if(isOK) { isOK = (StringMass > LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); }
871 ++cClusterInterrupt;
872 //... repeat procedure, if mass of cluster is too low to produce hadrons
873 //... ClusterMassCut = 0.15*GeV model parameter
874 }
875 while (isOK == false || cClusterInterrupt < ClusterLoopInterrupt);
876 /* Loop checking, 07.08.2015, A.Ribon */
877 return isOK;
878}
879
880//----------------------------------------------------------------------------------------
881
884 G4ParticleDefinition * & LeftHadron,
885 G4ParticleDefinition * & RightHadron)
886{
887 // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
888
889 G4double StringMass = string->Mass();
890 G4double StringMassSqr= sqr(StringMass);
891 G4ParticleDefinition * Di_Quark;
892 G4ParticleDefinition * Anti_Di_Quark;
893
894 if (string->GetLeftParton()->GetPDGEncoding() < 0)
895 {
896 Anti_Di_Quark =string->GetLeftParton();
897 Di_Quark=string->GetRightParton();
898 } else
899 {
900 Anti_Di_Quark =string->GetRightParton();
901 Di_Quark=string->GetLeftParton();
902 }
903
904 G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
905 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
906 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
907 G4int AbsIDdi_quark =std::abs(IDdi_quark);
908
909 G4int ADi_q1=AbsIDAnti_di_quark/1000;
910 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
911
912 G4int Di_q1=AbsIDdi_quark/1000;
913 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
914
915 NumberOf_FS=0;
916 for (G4int ProdQ=1; ProdQ < 6; ProdQ++)
917 {
918 G4int StateADiQ=0;
919 const G4int maxNumberOfLoops = 1000;
920 G4int loopCounter = 0;
921 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
922 {
924 -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
925
926 if (LeftHadron == NULL) continue;
927 G4double LeftHadronMass=LeftHadron->GetPDGMass();
928
929 G4int StateDiQ=0;
930 const G4int maxNumberOfInternalLoops = 1000;
931 G4int internalLoopCounter = 0;
932 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
933 {
935 +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
936
937 if (RightHadron == NULL) continue;
938 G4double RightHadronMass=RightHadron->GetPDGMass();
939
940 if (StringMass > LeftHadronMass + RightHadronMass)
941 {
942 if ( NumberOf_FS > 349 ) {
944 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
945 G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
946 "HAD_LUND_001", JustWarning, ed );
947 NumberOf_FS = 349;
948 }
949
950 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
951 sqr(RightHadronMass));
952 //FS_Psqr=1.;
953 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
954 BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
955 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
956 Prob_QQbar[ProdQ-1];
957
958 FS_LeftHadron[NumberOf_FS] = LeftHadron;
959 FS_RightHadron[NumberOf_FS]= RightHadron;
960 NumberOf_FS++;
961 } // End of if (StringMass > LeftHadronMass + RightHadronMass)
962
963 StateDiQ++;
964
965 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
966 ++internalLoopCounter < maxNumberOfInternalLoops );
967 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
968 return false;
969 }
970
971 StateADiQ++;
972 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
973 ++loopCounter < maxNumberOfLoops );
974 if ( loopCounter >= maxNumberOfLoops ) {
975 return false;
976 }
977 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
978
979 return true;
980}
981
982//----------------------------------------------------------------------------------------
983
985 G4ParticleDefinition * & LeftHadron,
986 G4ParticleDefinition * & RightHadron)
987{
988 G4double StringMass = string->Mass();
989 G4double StringMassSqr= sqr(StringMass);
990
991 G4ParticleDefinition * Di_Quark;
992 G4ParticleDefinition * Quark;
993
994 if (string->GetLeftParton()->GetParticleSubType()== "quark")
995 {
996 Quark =string->GetLeftParton();
997 Di_Quark=string->GetRightParton();
998 } else
999 {
1000 Quark =string->GetRightParton();
1001 Di_Quark=string->GetLeftParton();
1002 }
1003
1004 G4int IDquark =Quark->GetPDGEncoding();
1005 G4int AbsIDquark =std::abs(IDquark);
1006 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1007 G4int AbsIDdi_quark=std::abs(IDdi_quark);
1008 G4int Di_q1=AbsIDdi_quark/1000;
1009 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1010 G4int SignDiQ= 1;
1011 if (IDdi_quark < 0) SignDiQ=-1;
1012
1013 NumberOf_FS=0;
1014 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar
1015 { // (as last splitting, do not consider c-cbar and b-bbar cases)
1016 G4int SignQ;
1017 if (IDquark > 0)
1018 {
1019 SignQ=-1;
1020 if (IDquark == 2) SignQ= 1;
1021 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1022 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1023 if (IDquark == 4) SignQ= 1; // D+, D0, Ds+
1024 if (IDquark == 5) SignQ=-1; // B-, anti_B0, anti_Bs0
1025 } else
1026 {
1027 SignQ= 1;
1028 if (IDquark == -2) SignQ=-1;
1029 if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1030 if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1031 if (IDquark == -4) SignQ=-1; // D-, anti_D0, anti_Ds+
1032 if (IDquark == -5) SignQ= 1; // B+, B0, Bs0
1033 }
1034
1035 if (AbsIDquark == ProdQ) SignQ= 1;
1036
1037 G4int StateQ=0;
1038 const G4int maxNumberOfLoops = 1000;
1039 G4int loopCounter = 0;
1040 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1041 {
1043 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1044 if (LeftHadron == NULL) continue;
1045 G4double LeftHadronMass=LeftHadron->GetPDGMass();
1046
1047 G4int StateDiQ=0;
1048 const G4int maxNumberOfInternalLoops = 1000;
1049 G4int internalLoopCounter = 0;
1050 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1051 {
1052 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1053 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1054 if (RightHadron == NULL) continue;
1055 G4double RightHadronMass=RightHadron->GetPDGMass();
1056
1057 if (StringMass > LeftHadronMass + RightHadronMass)
1058 {
1059 if ( NumberOf_FS > 349 ) {
1061 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1062 G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1063 "HAD_LUND_002", JustWarning, ed );
1064 NumberOf_FS = 349;
1065 }
1066
1067 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1068 sqr(RightHadronMass));
1069 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1070 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1071 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1072 Prob_QQbar[ProdQ-1];
1073
1074 FS_LeftHadron[NumberOf_FS] = LeftHadron;
1075 FS_RightHadron[NumberOf_FS]= RightHadron;
1076
1077 NumberOf_FS++;
1078 } // End of if (StringMass > LeftHadronMass + RightHadronMass)
1079
1080 StateDiQ++;
1081
1082 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1083 ++internalLoopCounter < maxNumberOfInternalLoops );
1084 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1085 return false;
1086 }
1087
1088 StateQ++;
1089 } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1090 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1091
1092 if ( loopCounter >= maxNumberOfLoops ) {
1093 return false;
1094 }
1095 }
1096
1097 return true;
1098}
1099
1100//----------------------------------------------------------------------------------------
1101
1103 G4ParticleDefinition * & LeftHadron,
1104 G4ParticleDefinition * & RightHadron)
1105{
1106 G4double StringMass = string->Mass();
1107 G4double StringMassSqr= sqr(StringMass);
1108
1109 G4ParticleDefinition * Quark;
1110 G4ParticleDefinition * Anti_Quark;
1111
1112 if (string->GetLeftParton()->GetPDGEncoding()>0)
1113 {
1114 Quark =string->GetLeftParton();
1115 Anti_Quark=string->GetRightParton();
1116 } else
1117 {
1118 Quark =string->GetRightParton();
1119 Anti_Quark=string->GetLeftParton();
1120 }
1121
1122 G4int IDquark =Quark->GetPDGEncoding();
1123 G4int AbsIDquark =std::abs(IDquark);
1124 G4int QuarkCharge =Qcharge[IDquark-1];
1125
1126 G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1127 G4int AbsIDanti_quark =std::abs(IDanti_quark);
1128 G4int AntiQuarkCharge =-Qcharge[AbsIDanti_quark-1];
1129
1130 G4int LeftHadronCharge(0), RightHadronCharge(0);
1131
1132 //G4cout<<"Q Qbar "<<IDquark<<" "<<IDanti_quark<<G4endl;
1133
1134 NumberOf_FS=0;
1135 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar
1136 { // (as last splitting, do not consider c-cbar and b-bbar cases)
1137 LeftHadronCharge = QuarkCharge - Qcharge[ProdQ-1];
1138 G4int SignQ = LeftHadronCharge/3; if (SignQ == 0) SignQ = 1;
1139
1140 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 (d,sbar)
1141 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar (s,dbar)
1142 if ((IDquark == 4) && (ProdQ == 2)) SignQ= 1; // D0 (c,ubar)
1143 if ((IDquark == 5) && (ProdQ == 1)) SignQ=-1; // anti_B0 (b,dbar)
1144 if ((IDquark == 5) && (ProdQ == 3)) SignQ=-1; // anti_Bs0 (b,sbar)
1145
1146 RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1];
1147 G4int SignAQ = RightHadronCharge/3; if (SignAQ == 0) SignAQ = 1;
1148
1149 if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar (dbar,s)
1150 if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0 (sbar,d)
1151 if ((IDanti_quark ==-4) && (ProdQ == 2)) SignAQ=-1; // anti_D0 (cbar,u)
1152 if ((IDanti_quark ==-5) && (ProdQ == 1)) SignAQ= 1; // B0 (bbar,d)
1153 if ((IDanti_quark ==-5) && (ProdQ == 3)) SignAQ= 1; // Bs0 (bbar,s)
1154
1155 //G4cout<<"ProQ signs "<<ProdQ<<" "<<SignQ<<" "<<SignAQ<<G4endl;
1156
1157 G4int StateQ=0;
1158 const G4int maxNumberOfLoops = 1000;
1159 G4int loopCounter = 0;
1160 do
1161 {
1162 //G4cout<<"[AbsIDquark-1][ProdQ-1][StateQ "<<AbsIDquark-1<<" "
1163 //<<ProdQ-1<<" "<<StateQ<<" "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<G4endl;
1165 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1166 //G4cout<<"LeftHadron "<<LeftHadron<<G4endl;
1167 if (LeftHadron == NULL) { StateQ++; continue; }
1168 //G4cout<<"LeftHadron "<<LeftHadron->GetParticleName()<<G4endl;
1169 G4double LeftHadronMass=LeftHadron->GetPDGMass();
1170
1171 G4int StateAQ=0;
1172 const G4int maxNumberOfInternalLoops = 1000;
1173 G4int internalLoopCounter = 0;
1174 do
1175 {
1176 //G4cout<<" [AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<AbsIDanti_quark-1<<" "
1177 // <<ProdQ-1<<" "<<StateAQ<<" "<<SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<G4endl;
1178 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1179 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1180 //G4cout<<"RightHadron "<<RightHadron<<G4endl;
1181 if(RightHadron == NULL) { StateAQ++; continue; }
1182 //G4cout<<"RightHadron "<<RightHadron->GetParticleName()<<G4endl;
1183 G4double RightHadronMass=RightHadron->GetPDGMass();
1184
1185 if (StringMass > LeftHadronMass + RightHadronMass)
1186 {
1187 if ( NumberOf_FS > 349 ) {
1189 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1190 G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1191 "HAD_LUND_003", JustWarning, ed );
1192 NumberOf_FS = 349;
1193 }
1194
1195 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1196 sqr(RightHadronMass));
1197 //FS_Psqr=1.;
1198 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1199 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1200 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1201 Prob_QQbar[ProdQ-1];
1202
1203 if (string->GetLeftParton()->GetPDGEncoding()>0)
1204 {
1205 FS_LeftHadron[NumberOf_FS] = RightHadron;
1206 FS_RightHadron[NumberOf_FS]= LeftHadron;
1207 } else
1208 {
1209 FS_LeftHadron[NumberOf_FS] = LeftHadron;
1210 FS_RightHadron[NumberOf_FS]= RightHadron;
1211 }
1212
1213 NumberOf_FS++;
1214
1215 }
1216
1217 StateAQ++;
1218 //G4cout<<" StateAQ Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<StateAQ<<" "
1219 // <<Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<" "<<internalLoopCounter<<G4endl;
1220 } while ( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1221 ++internalLoopCounter < maxNumberOfInternalLoops );
1222 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1223 return false;
1224 }
1225
1226 StateQ++;
1227 //G4cout<<"StateQ Meson[AbsIDquark-1][ProdQ-1][StateQ] "<<StateQ<<" "
1228 // <<Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<loopCounter<<G4endl;
1229
1230 } while ( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1231 ++loopCounter < maxNumberOfLoops );
1232 if ( loopCounter >= maxNumberOfLoops ) {
1233 return false;
1234 }
1235 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
1236
1237 return true;
1238}
1239
1240//----------------------------------------------------------------------------------------------------------
1241
1243{
1244 if ( NumberOf_FS > 349 ) {
1246 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1247 G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1248 NumberOf_FS = 349;
1249 }
1250
1251 G4double SumWeights=0.;
1252
1253 for (G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1254
1255 G4double ksi=G4UniformRand();
1256 G4double Sum=0.;
1257 G4int indexPosition = 0;
1258
1259 for (G4int i=0; i<NumberOf_FS; i++)
1260 {
1261 Sum+=(FS_Weight[i]/SumWeights);
1262 indexPosition=i;
1263 if (Sum >= ksi) break;
1264 }
1265 return indexPosition;
1266}
1267
1268//----------------------------------------------------------------------------------------------------------
1269
1271 G4LorentzVector* AntiMom, G4double AntiMass,
1272 G4double InitialMass)
1273{
1274 // ------ Sampling of momenta of 2 last produced hadrons --------------------
1275 G4ThreeVector Pt;
1276 G4double MassMt, AntiMassMt;
1277 G4double AvailablePz, AvailablePz2;
1278 #ifdef debug_LUNDfragmentation
1279 G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
1280 G4cout<<"Init Mass "<<InitialMass<<" FirstM "<<Mass<<" SecondM "<<AntiMass<<" ProbIsotropy "<<G4endl;
1281 #endif
1282
1283 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
1284 sqr(2.*Mass*AntiMass);
1285 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
1286
1287 const G4int maxNumberOfLoops = 1000;
1288 G4double SigmaQTw = SigmaQT;
1289 if ( Mass > 930. || AntiMass > 930. ) {
1290 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1291 }
1292 if ( Mass < 930. && AntiMass < 930. ) {} // q-qbar string
1293 if ( ( Mass < 930. && AntiMass > 930. ) ||
1294 ( Mass > 930. && AntiMass < 930. ) ) { // q-di_q string
1295 //SigmaQT = -1.; // isotropical decay
1296 }
1297 if ( Mass > 930. && AntiMass > 930. ) { // qq-qqbar string
1298 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1299 }
1300
1301 G4int loopCounter = 0;
1302 do
1303 {
1304 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2();
1305 MassMt = std::sqrt( Mass * Mass + Pt2);
1306 AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1307 }
1308 while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops );
1309
1310 SigmaQT = SigmaQTw;
1311
1312 if ( loopCounter >= maxNumberOfLoops ) {
1313 AvailablePz2 = 0.0;
1314 }
1315
1316 AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) -
1317 4.*sqr(MassMt*AntiMassMt);
1318
1319 AvailablePz2 /=(4.*InitialMass*InitialMass);
1320 AvailablePz = std::sqrt(AvailablePz2);
1321
1322 G4double Px=Pt.getX();
1323 G4double Py=Pt.getY();
1324
1325 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1326 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2));
1327
1328 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1329 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2));
1330
1331 #ifdef debug_LUNDfragmentation
1332 G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl;
1333 G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ()
1334 <<" "<<AntiMom->getT()<<G4endl;
1335 #endif
1336}
1337
1338//------------------------------------------------------------------------
1339
1341{
1342 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1343 return lam;
1344}
1345
1346// --------------------------------------------------------------
1348{}
1349
G4double S(G4double temp)
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4double G4Log(G4double x)
Definition: G4Log.hh:226
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double fermi
Definition: G4SIunits.hh:83
static constexpr double GeV
Definition: G4SIunits.hh:203
static constexpr double MeV
Definition: G4SIunits.hh:200
static constexpr double pi
Definition: G4SIunits.hh:55
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double x() const
double mag2() const
double y() const
void setZ(double)
double getX() const
double getY() const
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
double theta() const
Hep3Vector boostVector() const
Hep3Vector vect() const
G4double GetTimeOfCreation() const
const G4ThreeVector & GetPosition() const
G4int GetDirection(void) const
G4Parton * GetRightParton(void) const
G4Parton * GetLeftParton(void) const
G4LorentzVector Get4Momentum() const
G4bool IsAFourQuarkString(void) const
G4ParticleDefinition * GetLeftParton(void) const
G4ParticleDefinition * GetRightParton(void) const
void LorentzRotate(const G4LorentzRotation &rotation)
G4LorentzRotation TransformToAlignedCms()
G4ParticleDefinition * GetDecayParton() const
G4int GetDecayDirection() const
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
static G4HadronicParameters * Instance()
G4double GetFormationTime() const
void SetPosition(const G4ThreeVector aPosition)
void Set4Momentum(const G4LorentzVector &a4Momentum)
const G4ThreeVector & GetPosition() const
const G4ParticleDefinition * GetDefinition() const
void SetFormationTime(G4double aFormationTime)
const G4LorentzVector & Get4Momentum() const
G4bool Quark_AntiQuark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
G4bool Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
G4bool Quark_Diquark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
virtual G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
G4double lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
G4bool Loop_toFragmentString(const G4ExcitedString &theStringInCMS, G4KineticTrackVector *&LeftVector, G4KineticTrackVector *&RightVector)
G4bool Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
virtual G4bool IsItFragmentable(const G4FragmentingString *const string)
const G4String & GetParticleName() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4int GetPDGcode() const
Definition: G4Parton.hh:127
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
G4ParticleDefinition * FindParticle(G4int Encoding)
G4KineticTrackVector * ProduceOneHadron(const G4ExcitedString *const theString)
virtual G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
G4ParticleDefinition * FS_RightHadron[350]
virtual void SetMassCut(G4double aValue)
void SetDiquarkSuppression(G4double aValue)
void SetStrangenessSuppression(G4double aValue)
G4ParticleDefinition * FS_LeftHadron[350]
void SetDiquarkBreakProbability(G4double aValue)
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
void SetStringTensionParameter(G4double aValue)
void SetMinimalStringMass(const G4FragmentingString *const string)
ush Pos
Definition: deflate.h:91
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4int sign(const T t)
float c_light
Definition: hepunit.py:256
T sqr(const T &x)
Definition: templates.hh:128