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

Last change on this file since 978 was 962, checked in by garnier, 15 years ago

update processes
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-ref-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.