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

Last change on this file since 1245 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
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 $ 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.