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

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 28.4 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4LundStringFragmentation.cc,v 1.13 2008/06/23 09:17:10 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02 $ 1.8
29//
30// -----------------------------------------------------------------------------
31// GEANT 4 class implementation file
32//
33// History: first implementation, Maxim Komogorov, 10-Jul-1998
34// -----------------------------------------------------------------------------
35#include "G4LundStringFragmentation.hh"
36#include "G4FragmentingString.hh"
37#include "G4DiQuarks.hh"
38#include "G4Quarks.hh"
39
40#include "Randomize.hh"
41
42// Class G4LundStringFragmentation
43//*************************************************************************************
44
45G4LundStringFragmentation::G4LundStringFragmentation()
46 {
47// ------ For estimation of a minimal string mass ---------------
48 Mass_of_light_quark =140.*MeV;
49 Mass_of_heavy_quark =500.*MeV;
50 Mass_of_string_junction=720.*MeV;
51// ------ An estimated minimal string mass ----------------------
52 MinimalStringMass = 0.;
53 MinimalStringMass2 = 0.;
54// ------ Minimal invariant mass used at a string fragmentation -
55 WminLUND = 0.7*GeV; // Uzhi 0.8 1.5
56// ------ Smooth parameter used at a string fragmentation for ---
57// ------ smearinr sharp mass cut-off ---------------------------
58 SmoothParam = 0.2;
59
60 SetStringTensionParameter(0.25); // Uzhi 20 June 08
61 }
62
63// --------------------------------------------------------------
64G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
65 {
66 }
67
68G4LundStringFragmentation::~G4LundStringFragmentation()
69 {
70 }
71
72//*************************************************************************************
73
74const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
75 {
76 throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
77 return *this;
78 }
79
80int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
81 {
82 return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
83 }
84
85int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
86 {
87 return memcmp(this, &right, sizeof(G4LundStringFragmentation));
88 }
89
90//--------------------------------------------------------------------------------------
91void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string) // Uzhi
92{
93/*
94G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
95G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<
96 string->GetRightParton()->GetPDGEncoding()<<G4endl;
97*/
98 G4double EstimatedMass=0.;
99 G4int Number_of_quarks=0;
100
101 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
102
103 if( Qleft > 1000)
104 {
105 Number_of_quarks+=2;
106 G4int q1=Qleft/1000;
107 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
108 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
109
110 G4int q2=(Qleft/100)%10;
111 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
112 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
113 EstimatedMass +=Mass_of_string_junction;
114 }
115 else
116 {
117 Number_of_quarks++;
118 if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
119 if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;}
120 }
121
122 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
123
124 if( Qright > 1000)
125 {
126 Number_of_quarks+=2;
127 G4int q1=Qright/1000;
128 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
129 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
130
131 G4int q2=(Qright/100)%10;
132 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
133 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
134 EstimatedMass +=Mass_of_string_junction;
135 }
136 else
137 {
138 Number_of_quarks++;
139 if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
140 if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;}
141 }
142
143 if(Number_of_quarks==2){EstimatedMass +=100.*MeV;}
144 if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
145 if(Number_of_quarks==4){EstimatedMass -=2.*Mass_of_string_junction;
146 if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
147 else {EstimatedMass+=100.*MeV;}
148 }
149 MinimalStringMass=EstimatedMass;
150 SetMinimalStringMass2(EstimatedMass);
151//G4cout<<"Out SetMinimalStringMass "<<MinimalStringMass<<G4endl;
152}
153
154//--------------------------------------------------------------------------------------
155void G4LundStringFragmentation::SetMinimalStringMass2(
156 const G4double aValue)
157{
158 MinimalStringMass2=aValue * aValue;
159}
160
161//--------------------------------------------------------------------------------------
162G4KineticTrackVector* G4LundStringFragmentation::FragmentString(
163 const G4ExcitedString& theString)
164{
165
166//G4cout<<"In FragmentString"<<G4endl;
167
168// Can no longer modify Parameters for Fragmentation.
169 PastInitPhase=true;
170
171// check if string has enough mass to fragment...
172
173 SetMassCut(160.*MeV); // For LightFragmentationTest it is required
174 // that no one pi-meson can be produced
175/*
176G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl;
177G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<
178theString.GetTimeOfCreation()/fermi<<G4endl;
179G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl;
180*/
181 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
182 if ( LeftVector != 0 ) {
183//G4cout<<"Return single hadron insted of string"<<G4endl;
184// Uzhi insert 6.05.08 start
185 if(LeftVector->size() == 1){
186 // One hadron is saved in the interaction
187 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
188 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
189
190/* // To set large formation time open *
191 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()+100.*fermi);
192 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
193 G4ThreeVector aPosition(theString.GetPosition().x(),
194 theString.GetPosition().y(),
195 theString.GetPosition().z()+100.*fermi);
196 LeftVector->operator[](0)->SetPosition(aPosition);
197*/
198//G4cout<<"Single hadron "<<LeftVector->operator[](0)->GetPosition()<<" "<<LeftVector->operator[](0)->GetFormationTime()<<G4endl;
199 } else { // 2 hadrons created from qq-qqbar are stored
200 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
201 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
202 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
203 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
204 }
205// Uzhi insert 6.05.08 end
206 return LeftVector;
207 }
208
209//--------------------- The string can fragment -------------------------------
210//--------------- At least two particles can be produced ----------------------
211
212 LeftVector =new G4KineticTrackVector;
213 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
214
215 G4ExcitedString *theStringInCMS=CPExcited(theString);
216 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
217
218 G4bool success=false, inner_sucess=true;
219 G4int attempt=0;
220 while ( !success && attempt++ < StringLoopInterrupt )
221 { // If the string fragmentation do not be happend, repeat the fragmentation---
222 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
223
224//G4cout<<"FragmentString cur M2 "<<std::sqrt(currentString->Mass2())<<G4endl;
225 // Cleaning up the previously produced hadrons ------------------------------
226 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
227 LeftVector->clear();
228
229 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
230 RightVector->clear();
231
232 // Main fragmentation loop until the string will not be able to fragment ----
233 inner_sucess=true; // set false on failure..
234 while (! StopFragmenting(currentString) )
235 { // Split current string into hadron + new string
236 G4FragmentingString *newString=0; // used as output from SplitUp...
237
238//G4cout<<"FragmentString to Splitup ===================================="<<G4endl;
239//G4cout<<"++++++++++++++++++++++++++ Enter num--------------------------"<<G4endl;
240//G4int Uzhi; G4cin>>Uzhi; // Uzhi
241
242 G4KineticTrack * Hadron=Splitup(currentString,newString);
243//G4cout<<" Hadron "<<Hadron<<G4endl;
244
245 if ( Hadron != 0 ) // Store the hadron
246 {
247 if ( currentString->GetDecayDirection() > 0 )
248 LeftVector->push_back(Hadron);
249 else
250 RightVector->push_back(Hadron);
251 delete currentString;
252 currentString=newString;
253 }
254 };
255 // Split remaining string into 2 final Hadrons ------------------------
256//G4cout<<"FragmentString to SplitLast if inner_sucess#0"<<inner_sucess<<G4endl;
257 if ( inner_sucess &&
258 SplitLast(currentString,LeftVector, RightVector) )
259 {
260 success=true;
261 }
262 delete currentString;
263 } // End of the loop in attemps to fragment the string
264
265 delete theStringInCMS;
266
267 if ( ! success )
268 {
269 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
270 LeftVector->clear();
271 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
272 delete RightVector;
273 return LeftVector;
274 }
275
276 // Join Left- and RightVector into LeftVector in correct order.
277 while(!RightVector->empty())
278 {
279 LeftVector->push_back(RightVector->back());
280 RightVector->erase(RightVector->end()-1);
281 }
282 delete RightVector;
283
284//G4cout<<"CalculateHadronTimePosition"<<G4endl;
285
286 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
287
288 G4LorentzRotation toObserverFrame(toCms.inverse());
289
290// LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
291// LeftVector->operator[](0)->SetPosition(theString.GetPosition());
292
293 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
294 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
295
296/* // For large formation time open *
297 G4double TimeOftheStringCreation=theString.GetTimeOfCreation()+100*fermi;
298 G4ThreeVector PositionOftheStringCreation(theString.GetPosition().x(),
299 theString.GetPosition().y(),
300 theString.GetPosition().z()+100*fermi);
301*/
302
303/*
304 if(theString.GetPosition().y() > 100.*fermi){
305// It is a projectile-like string -------------------------------------
306 G4double Zmin=theString.GetPosition().y()-1000.*fermi;
307 G4double Zmax=theString.GetPosition().z();
308 TimeOftheStringCreation=
309 (Zmax-Zmin)*theString.Get4Momentum().e()/theString.Get4Momentum().z();
310
311 G4ThreeVector aPosition(0.,0.,Zmax);
312 PositionOftheStringCreation=aPosition;
313 }
314*/
315 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
316 {
317 G4KineticTrack* Hadron = LeftVector->operator[](C1);
318 G4LorentzVector Momentum = Hadron->Get4Momentum();
319 Momentum = toObserverFrame*Momentum;
320 Hadron->Set4Momentum(Momentum);
321
322 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
323 Momentum = toObserverFrame*Coordinate;
324 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e());
325 G4ThreeVector aPosition(Momentum.vect());
326// Hadron->SetPosition(theString.GetPosition()+aPosition);
327 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
328//G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl;
329 };
330
331//G4cout<<"Out FragmentString"<<G4endl;
332 return LeftVector;
333
334}
335
336//----------------------------------------------------------------------------------
337G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
338{
339//G4cout<<"In IsFragmentable"<<G4endl;
340 SetMinimalStringMass(string); // Uzhi
341//G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
342 return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2(); // Uzhi
343
344// return sqr(FragmentationMass(string)+MassCut) < // Uzhi
345// string->Mass2(); // Uzhi
346}
347
348//----------------------------------------------------------------------------------------
349G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
350{
351//G4cout<<"StopFragmenting"<<G4endl;
352
353 SetMinimalStringMass(string);
354
355//G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
356
357 return (MinimalStringMass + WminLUND)*
358 (1 + SmoothParam * (1.-2*G4UniformRand())) >
359 string->Mass();
360}
361
362//----------------------------------------------------------------------------------------
363G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
364 G4KineticTrackVector * LeftVector,
365 G4KineticTrackVector * RightVector)
366{
367 //... perform last cluster decay
368/*
369G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl;
370G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
371G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
372*/
373 G4LorentzVector Str4Mom=string->Get4Momentum();
374//G4cout<<"String 4 momentum "<<Str4Mom<<G4endl;
375 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
376
377 G4double ResidualMass = string->Mass();
378
379 G4ParticleDefinition * LeftHadron, * RightHadron;
380
381 G4int cClusterInterrupt = 0;
382 do
383 {
384//G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl;
385
386 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
387 {
388 return false;
389 }
390 G4ParticleDefinition * quark = NULL;
391 string->SetLeftPartonStable(); // to query quark contents..
392
393 if (!string->FourQuarkString() )
394 {
395 // The string is q-qbar, or q-qq, or qbar-qqbar type
396 if (string->DecayIsQuark() && string->StableIsQuark() )
397 {
398 //... there are quarks on cluster ends
399 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
400 } else
401 {
402 //... there is a Diquark on one of the cluster ends
403 G4int IsParticle;
404
405 if ( string->StableIsQuark() )
406 {
407 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
408 } else
409 {
410 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
411 }
412
413 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
414 quark = QuarkPair.second;
415
416 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
417 }
418
419 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
420 } else
421 {
422 // The string is qq-qqbar type. Diquarks are on the string ends
423 G4int LiftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
424 G4int LiftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
425
426 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
427 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
428
429 if(G4UniformRand()<0.5)
430 {
431 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
432 FindParticle(RightQuark1));
433 RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
434 FindParticle(RightQuark2));
435 } else
436 {
437 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
438 FindParticle(RightQuark2));
439 RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
440 FindParticle(RightQuark1));
441 }
442 }
443/*
444G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl;
445G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl;
446G4cout<<"Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl;
447*/
448 //... repeat procedure, if mass of cluster is too low to produce hadrons
449 //... ClusterMassCut = 0.15*GeV model parameter
450
451 }
452 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());// UzhiVOVA
453
454 //... compute hadron momenta and energies
455
456 G4LorentzVector LeftMom, RightMom;
457 G4ThreeVector Pos;
458
459 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
460 &RightMom, RightHadron->GetPDGMass(),
461 ResidualMass);
462
463 LeftMom.boost(ClusterVel);
464 RightMom.boost(ClusterVel);
465
466 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
467 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
468
469//G4cout<<"Out SplitLast "<<G4endl;
470 return true;
471
472}
473
474//----------------------------------------------------------------------------------------------------------
475
476void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
477 {
478// ------ Sampling of momenta of 2 last produced hadrons --------------------
479 G4ThreeVector Pt;
480 G4double MassMt2, AntiMassMt2;
481 G4double AvailablePz, AvailablePz2;
482
483//G4cout<<"Sample4Momentum "<<G4endl;
484//G4cout<<"Sample4Momentum Mass"<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl;
485
486 if(Mass > 930. || AntiMass > 930.) // If there is a baryon
487 {
488 // ----------------- Isotripic decay ------------------------------------
489 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
490 sqr(2.*Mass*AntiMass);
491 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
492
493 //... sample unit vector
494 G4double pz = 1. - 2.*G4UniformRand();
495 G4double st = std::sqrt(1. - pz * pz)*Pabs;
496 G4double phi = 2.*pi*G4UniformRand();
497 G4double px = st*std::cos(phi);
498 G4double py = st*std::sin(phi);
499 pz *= Pabs;
500
501 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
502 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
503
504 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
505 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
506 }
507 else
508 {
509 do
510 {
511 Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2();
512
513//G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl;
514
515 MassMt2 = Mass * Mass + Pt2;
516 AntiMassMt2= AntiMass * AntiMass + Pt2;
517
518//G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl;
519
520 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
521 4.*MassMt2*AntiMassMt2;
522 }
523 while(AvailablePz2 < 0.);
524
525 AvailablePz2 /=(4.*InitialMass*InitialMass);
526 AvailablePz = std::sqrt(AvailablePz2);
527
528//G4cout<<"AvailablePz "<<AvailablePz<<G4endl;
529
530 G4double Px=Pt.getX();
531 G4double Py=Pt.getY();
532
533//if(Mass > AntiMass){AvailablePz=-AvailablePz;} // May30 // Uzhi
534
535 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
536 Mom->setE(std::sqrt(MassMt2+AvailablePz2));
537
538//G4cout<<" 1 part "<<Px<<" "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl;
539
540 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
541 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
542
543//G4cout<<" 2 part "<<-Px<<" "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl;
544 }
545
546//G4cout<<"Out Sample4Momentum "<<G4endl;
547 }
548
549//-----------------------------------------------------------------------------
550
551G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
552 G4FragmentingString * string, G4FragmentingString * newString)
553{
554/*
555G4cout<<"SplitEandP "<<G4endl;
556G4cout<<"SplitEandP string mass "<<string->Mass()<<G4endl;
557G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "
558 <<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
559G4cout<<G4endl;
560G4cout<<newString->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
561G4cout<<newString->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
562*/
563 G4LorentzVector String4Momentum=string->Get4Momentum();
564 G4double StringMT2=string->Get4Momentum().mt2();
565//G4cout<<"StringMt2 "<<StringMT2<<G4endl;
566
567 G4double HadronMass = pHadron->GetPDGMass();
568//G4cout<<"Hadron mass "<<HadronMass<<" "<<pHadron->GetParticleName()<<G4endl;
569
570 SetMinimalStringMass(newString);
571 String4Momentum.setPz(0.);
572 G4ThreeVector StringPt=String4Momentum.vect();
573//G4cout<<"StringPt "<<StringPt<<G4endl<<G4endl;
574//G4cout<<"Min string mass "<<MinimalStringMass<<G4endl;
575
576// calculate and assign hadron transverse momentum component HadronPx and HadronPy
577 G4ThreeVector thePt;
578 thePt=SampleQuarkPt();
579
580 G4ThreeVector HadronPt = thePt +string->DecayPt();
581 HadronPt.setZ(0);
582//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
583
584 G4ThreeVector RemSysPt = StringPt - HadronPt;
585//G4cout<<"RemSys Pt"<<RemSysPt<<G4endl;
586
587 //... sample z to define hadron longitudinal momentum and energy
588 //... but first check the available phase space
589
590 G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
591 G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
592
593//G4cout<<"Mt h res str "<<std::sqrt(HadronMassT2)<<" "<<std::sqrt(ResidualMassT2)<<" srt mass"<<StringMT2<<G4endl;
594
595 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
596 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
597
598//G4cout<<"Pz**2 "<<Pz2<<G4endl;
599
600 if(Pz2 < 0 ) {return 0;} // have to start all over!
601
602 //... then compute allowed z region z_min <= z <= z_max
603
604 G4double Pz = std::sqrt(Pz2);
605 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
606 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
607
608//G4cout<<"Zmin max "<<zMin<<" "<<zMax<<G4endl; // Uzhi
609
610 if (zMin >= zMax) return 0; // have to start all over!
611
612 G4double z = GetLightConeZ(zMin, zMax,
613 string->GetDecayParton()->GetPDGEncoding(), pHadron,
614 HadronPt.x(), HadronPt.y());
615
616 //... now compute hadron longitudinal momentum and energy
617 // longitudinal hadron momentum component HadronPz
618
619 HadronPt.setZ(0.5* string->GetDecayDirection() *
620 (z * string->LightConeDecay() -
621 HadronMassT2/(z * string->LightConeDecay())));
622
623 G4double HadronE = 0.5* (z * string->LightConeDecay() +
624 HadronMassT2/(z * string->LightConeDecay()));
625
626 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
627
628//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
629//G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<HadronE<<G4endl;
630
631 return a4Momentum;
632}
633
634//-----------------------------------------------------------------------------------------
635G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
636 G4int PDGEncodingOfDecayParton,
637 G4ParticleDefinition* pHadron,
638 G4double Px, G4double Py)
639{
640 G4double alund;
641
642// If blund get restored, you MUST adapt the calculation of zOfMaxyf.
643// const G4double blund = 1;
644
645 G4double z, yf;
646 G4double Mass = pHadron->GetPDGMass();
647// G4int HadronEncoding=pHadron->GetPDGEncoding();
648
649 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
650
651 if(std::abs(PDGEncodingOfDecayParton) < 1000)
652 {
653 // ---------------- Quark fragmentation ----------------------
654 alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered
655 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
656 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
657
658 do
659 {
660 z = zmin + G4UniformRand()*(zmax-zmin);
661// yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
662 yf = (1-z)/z * std::exp(-alund*Mt2/z);
663 }
664 while (G4UniformRand()*maxYf > yf);
665 }
666 else
667 {
668 // ---------------- Di-quark fragmentation ----------------------
669//G4cout<<"Di-quark"<<G4endl; // Vova
670 alund=0.7/GeV/GeV; // 0.7 2.0
671 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
672 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
673 do
674 {
675 z = zmin + G4UniformRand()*(zmax-zmin);
676// yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
677 yf = (1-z)/z * std::exp(-alund*Mt2/z);
678 }
679 while (G4UniformRand()*maxYf > yf);
680 };
681
682//G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl;
683 return z;
684 }
Note: See TracBrowser for help on using the repository browser.