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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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