| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4LundStringFragmentation.cc,v 1.5 2006/06/29 20:55:03 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: first implementation, Maxim K 33 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 34 // -----------------------------------------------------------------------------
33 #include "G4LundStringFragmentation.hh" 35 #include "G4LundStringFragmentation.hh"
34 #include "G4PhysicalConstants.hh" <<
35 #include "G4SystemOfUnits.hh" <<
36 #include "Randomize.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 37
46 // Class G4LundStringFragmentation 38 // Class G4LundStringFragmentation
47 //******************************************** << 39 //****************************************************************************************
48 40
49 G4LundStringFragmentation::G4LundStringFragmen 41 G4LundStringFragmentation::G4LundStringFragmentation()
50 : G4VLongitudinalStringDecay("LundStringFrag << 42 {
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 } 43 }
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 44
630 return a4Momentum; << 45 // G4LundStringFragmentation::G4LundStringFragmentation(G4double sigmaPt)
631 } << 46 // : G4VLongitudinalStringDecay(sigmaPt)
>> 47 // {
>> 48 // }
632 49
633 //-------------------------------------------- << 50 G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
634 << 51 {
635 G4double G4LundStringFragmentation::GetLightCo << 52 }
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 53
1327 SigmaQT = SigmaQTw; <<
1328 54
1329 AvailablePz2= sqr(InitialMass*InitialMass - << 55 G4LundStringFragmentation::~G4LundStringFragmentation()
1330 4.*sqr(MassMt*AntiMassMt); << 56 {
>> 57 }
1331 58
1332 AvailablePz2 /=(4.*InitialMass*InitialMass) << 59 //****************************************************************************************
1333 AvailablePz = std::sqrt(AvailablePz2); <<
1334 60
1335 G4double Px=Pt.getX(); << 61 const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
1336 G4double Py=Pt.getY(); << 62 {
>> 63 throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
>> 64 return *this;
>> 65 }
1337 66
1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( << 67 int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz << 68 {
>> 69 return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
>> 70 }
1340 71
1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A << 72 int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av << 73 {
>> 74 return memcmp(this, &right, sizeof(G4LundStringFragmentation));
>> 75 }
1343 76
1344 #ifdef debug_LUNDfragmentation << 77 //****************************************************************************************
1345 G4cout<<"Fmass Mom "<<Mom->getX()<<" <<
1346 G4cout<<"Smass Mom "<<AntiMom->getX() <<
1347 <<" "<<AntiMom->getT()<<G4endl; <<
1348 #endif <<
1349 } <<
1350 78
1351 //------------------------------------------- << 79 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int , G4ParticleDefinition* pHadron, G4double Px, G4double Py)
>> 80 {
>> 81 const G4double alund = 0.7/GeV/GeV;
1352 82
1353 G4double G4LundStringFragmentation::lambda(G4 << 83 // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
1354 { << 84 // const G4double blund = 1;
1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 <<
1356 return lam; <<
1357 } <<
1358 85
1359 // ------------------------------------------ << 86 G4double z, yf;
1360 G4LundStringFragmentation::~G4LundStringFragm << 87 G4double Mass = pHadron->GetPDGMass();
1361 {} << 88
>> 89 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
>> 90 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
>> 91 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
>> 92 do
>> 93 {
>> 94 z = zmin + G4UniformRand()*(zmax-zmin);
>> 95 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
>> 96 yf = (1-z)/z * std::exp(-alund*Mt2/z);
>> 97 }
>> 98 while (G4UniformRand()*maxYf > yf);
>> 99 return z;
>> 100 }
1362 101
>> 102 //****************************************************************************************
1363 103