Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc

    r819 r962  
    2525//
    2626//
    27 // $Id: G4LundStringFragmentation.cc,v 1.7 2007/04/24 14:55:23 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     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
    2929//
    3030// -----------------------------------------------------------------------------
     
    4141
    4242// Class G4LundStringFragmentation
    43 //****************************************************************************************
     43//*************************************************************************************
    4444
    4545G4LundStringFragmentation::G4LundStringFragmentation()
    4646   {
    47    MinimalStringMass  = 0.;             // Uzhi
    48    MinimalStringMass2 = 0.;             // Uzhi
    49    WminLUND = 1.*GeV;                   // Uzhi
    50    SmoothParam  = 0.2;                  // Uzhi
    51 
     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
    5261   }
    5362
    54 // G4LundStringFragmentation::G4LundStringFragmentation(G4double sigmaPt)
    55 // : G4VLongitudinalStringDecay(sigmaPt)
    56 //    {
    57 //    }
    58 
     63// --------------------------------------------------------------
    5964G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
    6065   {
    6166   }
    6267
    63 
    6468G4LundStringFragmentation::~G4LundStringFragmentation()
    6569   {
    6670   }
    6771
    68 //****************************************************************************************
     72//*************************************************************************************
    6973
    7074const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
     
    8488   }
    8589
    86 //****************************************************************************************
    87 //----------------------------------------------------------------------------------------------------------
    88 
    89 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
     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)
    90164{
    91165
     
    96170       
    97171//      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*/
    98181        G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
    99182        if ( LeftVector != 0 ) {
    100 //G4cout<<"Return single hadron from string"<<G4endl;
    101                                 return LeftVector;}
    102        
    103         LeftVector = new G4KineticTrackVector;
     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;
    104213        G4KineticTrackVector * RightVector=new G4KineticTrackVector;
    105214
    106 // this should work but its only a semi deep copy. %GF  G4ExcitedString theStringInCMS(theString);
    107215        G4ExcitedString *theStringInCMS=CPExcited(theString);
    108216        G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
     
    111219        G4int attempt=0;
    112220        while ( !success && attempt++ < StringLoopInterrupt )
    113         {
     221        { // If the string fragmentation do not be happend, repeat the fragmentation---
    114222                G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
    115223
    116 //G4cout<<"Main FragmentString cur M2  "<<std::sqrt(currentString->Mass2())<<G4endl;
    117 
    118                 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
     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());
    119227                LeftVector->clear();
     228
    120229                std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
    121230                RightVector->clear();
    122231               
     232          // Main fragmentation loop until the string will not be able to fragment ----
    123233                inner_sucess=true;  // set false on failure..
    124234                while (! StopFragmenting(currentString) )
    125235                {  // Split current string into hadron + new string
    126 
    127 //                      G4FragmentingString *PreviousString=currentString;   // Uzhi
    128 
    129236                        G4FragmentingString *newString=0;  // used as output from SplitUp...
    130237
    131238//G4cout<<"FragmentString to Splitup ===================================="<<G4endl;
     239//G4cout<<"++++++++++++++++++++++++++ Enter num--------------------------"<<G4endl;
    132240//G4int Uzhi; G4cin>>Uzhi;                         // Uzhi
     241
    133242                        G4KineticTrack * Hadron=Splitup(currentString,newString);
    134243//G4cout<<" Hadron "<<Hadron<<G4endl;
    135244
    136 //                      if ( Hadron != 0 && IsFragmentable(newString))       // Uzhi
    137                         if ( Hadron != 0 )                                   // Uzhi
     245                        if ( Hadron != 0 )  // Store the hadron                               
    138246                        {
    139247                           if ( currentString->GetDecayDirection() > 0 )
     
    143251                           delete currentString;
    144252                           currentString=newString;
    145                         } /* else {                                // Uzhi
    146                          // abandon ... start from the beginning
    147                            if (newString) delete newString;                  // ??? Uzhi local?
    148                            if (Hadron)    delete Hadron;
    149 //                           currentString = PreviousString;                 // Uzhi
    150                            inner_sucess=false;
    151                            break;
    152                         } */                                      // Uzhi
    153 //                        delete PreviousString;                             // ??? Uzhi local?
     253                        }
    154254                };
    155                 // Split current string into 2 final Hadrons
     255        // Split remaining string into 2 final Hadrons ------------------------
    156256//G4cout<<"FragmentString to SplitLast if inner_sucess#0"<<inner_sucess<<G4endl;
    157                 if ( inner_sucess &&                                         // Uzhi
     257                if ( inner_sucess &&                                       
    158258                     SplitLast(currentString,LeftVector, RightVector) )
    159259                {
     
    161261                }
    162262                delete currentString;
    163         }
     263        }  // End of the loop in attemps to fragment the string
    164264       
    165265        delete theStringInCMS;
     
    188288        G4LorentzRotation toObserverFrame(toCms.inverse());
    189289
     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*/
    190315        for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
    191316        {
     
    194319           Momentum = toObserverFrame*Momentum;
    195320           Hadron->Set4Momentum(Momentum);
     321
    196322           G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
    197323           Momentum = toObserverFrame*Coordinate;
    198            Hadron->SetFormationTime(Momentum.e());
     324           Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e());
    199325           G4ThreeVector aPosition(Momentum.vect());
    200            Hadron->SetPosition(theString.GetPosition()+aPosition);
    201         }
     326//         Hadron->SetPosition(theString.GetPosition()+aPosition);
     327           Hadron->SetPosition(PositionOftheStringCreation+aPosition);
     328//G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl;
     329        };
    202330
    203331//G4cout<<"Out FragmentString"<<G4endl;
    204332        return LeftVector;
    205333               
    206 
    207 
    208334}
    209335
     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
    210474//----------------------------------------------------------------------------------------------------------
    211475
    212 //G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,           // Uzhi
    213 //                                                  G4int ,  G4ParticleDefinition* pHadron, // Uzhi
    214 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
    215                                                   G4int,  G4ParticleDefinition* pHadron, // Uzhi
    216 G4double Px, G4double Py)
     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
    217487    {
    218     const G4double  alund = 0.7/GeV/GeV;
    219 
    220 //    If blund get restored, you MUST adapt the calculation of zOfMaxyf.
    221 //    const G4double  blund = 1;
    222 
    223     G4double z, yf;
    224     G4double Mass = pHadron->GetPDGMass();
     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;
    225500   
    226     G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
    227     G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
    228     G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
    229 
    230 //    G4double N=1.;                                                 // Uzhi
    231 //    G4double OverN=1./N;                                           // Uzhi
    232 //    G4double ZminN=std::pow(zmin,N);                               // Uzhi
    233 //    G4double ZmaxN=std::pow(zmax,N);                               // Uzhi
    234 //    G4double Brac=ZmaxN-ZminN;                                     // Uzhi
    235 
    236 //G4cout<<" ZminN ZmaxN Brac Code "<<ZminN<<" "<< ZmaxN<<" "<<Brac<<" "<<PartonEncoding<<G4endl;
    237 
    238 //    if(std::abs(PartonEncoding) < 1000)                            // Uzhi
    239       {                                                            // Uzhi q or q-bar
    240 //G4cout<<" quark "<<G4endl; // Vova
    241        do                                                          // Uzhi
    242          {
    243           z = zmin + G4UniformRand()*(zmax-zmin);
    244 //        yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
    245           yf = (1-z)/z * std::exp(-alund*Mt2/z);
    246          }
    247        while (G4UniformRand()*maxYf > yf);
    248       }                                                            // Uzhi
    249 //    else                                                           // Uzhi
    250 //      {                                                            // Uzhi qq or qq-bar
    251 // //G4cout<<"Di-quark"<<G4endl; // Vova
    252 //       z = std::pow(Brac * G4UniformRand() + ZminN, OverN);        // Uzhi
    253 //      };                                                           // Uzhi
    254 //
    255 //G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl; // Vova
    256     return z;
     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;
    257544    }
    258 //-----------------------------------------------------------------------------------------
     545
     546//G4cout<<"Out Sample4Momentum "<<G4endl;
     547  }
     548
     549//-----------------------------------------------------------------------------
    259550
    260551G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
    261         G4FragmentingString * string)
    262 {       
     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
    263567       G4double HadronMass = pHadron->GetPDGMass();
    264        SetMinimalStringMass(string);                                             // Uzhi
    265        G4double  StringMass2 = string->Mass2();                                  // Uzhi
    266 
    267 //G4cout<<"SplitEandP string mass "<<string->Mass()<<" Hadron mass "<<HadronMass<<pHadron->GetParticleName()<<G4endl;   // Uzhi
    268 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
    269 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
    270 //G4cout<<" Min string mass "<<MinimalStringMass<<G4endl;
    271 
    272        // calculate and assign hadron transverse momentum component HadronPx andHadronPy
     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
    273577       G4ThreeVector thePt;
    274578       thePt=SampleQuarkPt();
     
    276580       G4ThreeVector HadronPt = thePt +string->DecayPt();
    277581       HadronPt.setZ(0);
     582//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
     583
     584       G4ThreeVector RemSysPt = StringPt - HadronPt;
     585//G4cout<<"RemSys Pt"<<RemSysPt<<G4endl;
     586
    278587       //...  sample z to define hadron longitudinal momentum and energy
    279588       //... but first check the available phase space
    280589
    281 //       G4double DecayQuarkMass2  = sqr(string->GetDecayParton()->GetPDGMass());
    282 
    283 //G4cout<<" QuarkMass "<<string->GetDecayParton()->GetPDGMass()<<G4endl;              // Uzhi
    284 
    285        G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
    286 //       G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi
    287        G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi
    288 
    289 //G4cout<<" Mt h res str "<<std::sqrt(HadronMass2T)<<" "<<std::sqrt(ResidualMass2T)<<" srt mass"<<string->Mass()<<G4endl;
    290 
    291 //       if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )      // Uzhi
    292 
    293         G4double Pz2 = (sqr(StringMass2 - HadronMass2T - ResidualMass2T) -                  // Uzhi
    294                                         4*HadronMass2T * ResidualMass2T)/4./StringMass2; // Uzhi
    295 
    296 //G4cout<<" Pz**2 "<<Pz2<<G4endl;
    297 
    298         if(Pz2 < 0 ) {return 0;}          // have to start all over!                // Uzhi
     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!           
    299601
    300602       //... then compute allowed z region  z_min <= z <= z_max
    301603 
    302        G4double Pz = std::sqrt(Pz2);                                               // Uzhi
    303        G4double zMin = (std::sqrt(HadronMass2T+Pz2) - Pz)/std::sqrt(StringMass2);  // Uzhi
    304        G4double zMax = (std::sqrt(HadronMass2T+Pz2) + Pz)/std::sqrt(StringMass2);  // Uzhi
    305 
    306 //G4cout<<" Zmin max "<<zMin<<" "<<zMax<<G4endl;               // Uzhi
    307 
    308 //       G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());                   // Uzhi
     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
    309610       if (zMin >= zMax) return 0;              // have to start all over!
    310611       
     
    318619        HadronPt.setZ(0.5* string->GetDecayDirection() *
    319620                        (z * string->LightConeDecay() -
    320                          HadronMass2T/(z * string->LightConeDecay())));
     621                         HadronMassT2/(z * string->LightConeDecay())));
    321622
    322623        G4double HadronE  = 0.5* (z * string->LightConeDecay() +
    323                                   HadronMass2T/(z * string->LightConeDecay()));
     624                                  HadronMassT2/(z * string->LightConeDecay()));
    324625
    325626       G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
    326627
    327 //G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<0.5* (z * string->LightConeDecay() + HadronMass2T/(z * string->LightConeDecay()))<<G4endl;
     628//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
     629//G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<HadronE<<G4endl;
    328630
    329631       return a4Momentum;
    330632}
    331633
    332 
    333634//-----------------------------------------------------------------------------------------
    334 
    335 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
    336                                              G4KineticTrackVector * LeftVector,
    337                                              G4KineticTrackVector * RightVector)
     635G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
     636                                                  G4int PDGEncodingOfDecayParton,
     637                                                  G4ParticleDefinition* pHadron,
     638                                                  G4double Px, G4double Py)
    338639{
    339     //... perform last cluster decay
    340 //G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl;
    341 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
    342 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
    343 
    344     G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
    345 
    346     G4double ResidualMass   = string->Mass();
    347 //    G4double ClusterMassCut = ClusterMass;
    348 
    349     G4int cClusterInterrupt = 0;
    350 
    351     G4ParticleDefinition * LeftHadron, * RightHadron;
    352     do
    353     {
    354 //G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl;
    355 
    356         if (cClusterInterrupt++ >= ClusterLoopInterrupt)
    357         {
    358           return false;
    359         }
    360         G4ParticleDefinition * quark = NULL;
    361         string->SetLeftPartonStable(); // to query quark contents..
    362 
    363         if (string->DecayIsQuark() && string->StableIsQuark() )
    364         {
    365            //... there are quarks on cluster ends
    366                 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
    367         } else {
    368            //... there is a Diquark on cluster ends
    369                 G4int IsParticle;
    370 
    371                 if ( string->StableIsQuark() ) {
    372                   IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
    373                 } else {
    374                   IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
    375                 }
    376 
    377                 pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
    378                 quark = QuarkPair.second;
    379 
    380                 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
    381         }
    382 
    383         RightHadron = hadronizer->Build(string->GetRightParton(), quark);
    384 
    385 //G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl;
    386 //G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl;
    387 //G4cout<<" Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl;
    388 
    389        //... repeat procedure, if mass of cluster is too low to produce hadrons
    390        //... ClusterMassCut = 0.15*GeV model parameter
    391 //      if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}        // Uzhi
    392 //      else {ClusterMassCut = ClusterMass;}                                       // Uzhi
    393 
    394     }
    395     while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());         // Uzhi   VOVA
    396 //    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()  + ClusterMassCut); // Uzhi
    397     //... compute hadron momenta and energies   
    398 
    399     G4LorentzVector  LeftMom, RightMom;
    400     G4ThreeVector    Pos;
    401 
    402 //G4cout<<"Sample4Momentum"<<G4endl;
    403 
    404     Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
    405 
    406     LeftMom.boost(ClusterVel);
    407     RightMom.boost(ClusterVel);
    408 
    409     LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
    410     RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
    411 
    412     return true;
    413 
    414 }
    415 //----------------------------------------------------------------------------------------------------------
    416 
    417 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
    418 {
    419 //G4cout<<"In IsFragmentable"<<G4endl;
    420   SetMinimalStringMass(string);                                                            // Uzhi
    421 //G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
    422   return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();                // Uzhi
    423 
    424 //      return sqr(FragmentationMass(string)+MassCut) <                                    // Uzhi
    425 //                      string->Mass2();                                                   // Uzhi
    426 }
    427 
    428 //----------------------------------------------------------------------------------------------------------
    429 
    430 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
    431 {
    432 //G4cout<<"StopFragmenting"<<G4endl;
    433 
    434   SetMinimalStringMass(string);                                                            // Uzhi
    435 //G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
    436   return sqr((MinimalStringMass + WminLUND)*(1 + SmoothParam * (1.-2*G4UniformRand()))) >        // Uzhi
    437                    string->Get4Momentum().mag2();                                          // Uzhi
    438 //       sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >          // Uzhi
    439 //       string->Get4Momentum().mag2();                                                    // Uzhi
    440 }
    441 
    442 //----------------------------------------------------------------------------------------------------------
    443 
    444 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
    445     {
    446      G4ThreeVector Pt;                                                      // Uzhi
    447      G4double MassMt2, AntiMassMt2;                                         // Uzhi
    448      G4double AvailablePz, AvailablePz2;                                    // Uzhi
    449 
    450 //G4cout<<" Smpl4Mom "<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl;
    451                                                                             // Uzhi
    452     do                                                                      // Uzhi
    453       {                                                                     // Uzhi
    454        Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2();              // Uzhi
    455 
    456 //G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl;
    457 
    458        MassMt2    =     Mass *     Mass + Pt2;                              // Uzhi
    459        AntiMassMt2= AntiMass * AntiMass + Pt2;                              // Uzhi
    460 
    461 //G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl;
    462 
    463        AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
    464                      4.*MassMt2*AntiMassMt2;                                // Uzhi
    465       }                                                                     // Uzhi
    466     while(AvailablePz2 < 0.);                                               // Uzhi
    467                                                                             // Uzhi
    468     AvailablePz2 /=(4.*InitialMass*InitialMass);                            // Uzhi
    469                                                                             // Uzhi
    470     AvailablePz = std::sqrt(AvailablePz2);                                  // Uzhi
    471 
    472 //G4cout<<"AvailablePz "<<AvailablePz<<G4endl;
    473  
    474 
    475     G4double Px=Pt.getX();                                                  // Uzhi
    476     G4double Py=Pt.getY();                                                  // Uzhi
    477                                                                             // Uzhi
    478     Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);                // Uzhi
    479     Mom->setE(std::sqrt(MassMt2+AvailablePz2));                             // Uzhi
    480 
    481 //G4cout<<" 1 part "<<Px<<"  "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl;
    482 
    483                                                                             // Uzhi
    484     AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz); // Uzhi
    485     AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));                    // Uzhi
    486 
    487 //G4cout<<" 2 part "<<-Px<<"  "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl;
    488 
    489 // Maybe it must be inversed!                                               // Uzhi
    490 /*                                                                          // Uzhi
    491     G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
    492     G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
    493 
    494     //... sample unit vector       
    495     G4double pz = 1. - 2.*G4UniformRand(); 
    496     G4double st     = std::sqrt(1. - pz * pz)*Pabs;
    497     G4double phi    = 2.*pi*G4UniformRand();
    498     G4double px = st*std::cos(phi);
    499     G4double py = st*std::sin(phi);
    500     pz *= Pabs;
     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();
    501648   
    502     Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
    503     Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
    504 
    505     AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
    506     AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
    507 */                                                                           // Uzhi
     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;
    508684    }
    509 
    510 
    511 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString  * const string)  // Uzhi
    512 {
    513 //G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
    514 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<G4endl;
    515 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<G4endl;
    516 
    517   G4double EstimatedMass=0.750* GeV;  // 2*m_q
    518 
    519   G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
    520 
    521   if( Qleft > 1000)
    522     {
    523      G4int q1=Qleft/1000;
    524      if( q1 < 3) {EstimatedMass += 0.325* GeV;}
    525      if( q1 > 2) {EstimatedMass += 0.500* GeV;}     
    526 
    527      G4int q2=(Qleft/100)%10;
    528      if( q2 < 3) {EstimatedMass += 0.325* GeV;}
    529      if( q2 > 2) {EstimatedMass += 0.500* GeV;}
    530     }
    531   else
    532     {
    533      if( Qleft < 3) {EstimatedMass += 0.325* GeV;}
    534      if( Qleft > 2) {EstimatedMass += 0.500* GeV;}
    535     }
    536 
    537   G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
    538 
    539   if( Qright > 1000)
    540     {
    541      G4int q1=Qright/1000;
    542      if( q1 < 3) {EstimatedMass += 0.325* GeV;}
    543      if( q1 > 2) {EstimatedMass += 0.500* GeV;}     
    544 
    545      G4int q2=(Qright/100)%10;
    546      if( q2 < 3) {EstimatedMass += 0.325* GeV;}
    547      if( q2 > 2) {EstimatedMass += 0.500* GeV;}
    548     }
    549   else
    550     {
    551      if( Qright < 3) {EstimatedMass += 0.325* GeV;}
    552      if( Qright > 2) {EstimatedMass += 0.500* GeV;}
    553     }
    554 
    555   MinimalStringMass=EstimatedMass;
    556   SetMinimalStringMass2(EstimatedMass);
    557 
    558 /*
    559   Pcreate build=&G4HadronBuilder::BuildLowSpin;
    560 
    561   G4ParticleDefinition *Hadron1, *Hadron2=0;
    562 
    563   G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
    564 
    565   if (string->GetLeftParton()->GetParticleSubType() == "quark") iflc = -iflc;
    566 
    567   if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
    568 
    569   // 1/2 baryon (anti-baryon) and scalar meson     (QQ-q or QbarQbar-Qbar),
    570   // or 2 scalar mesons                            (Q-Qbar),
    571   // or 2 1/2 baryons (anti-baryons) will be built (QQ-QbarQbar)
    572 
    573 //G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
    574 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<FindParticle(iflc)->GetPDGEncoding()<<G4endl;
    575 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<FindParticle(-iflc)->GetPDGEncoding()<<G4endl;
    576 
    577   Hadron1 = (hadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));
    578   Hadron2 =(hadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));
    579   MinimalStringMass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
    580 
    581 //G4cout<<(Hadron1)->GetPDGEncoding()<<" "<<(Hadron2)->GetPDGEncoding()<<G4endl;
    582 //G4cout<<"Out SetMinMass "<<MinimalStringMass<<G4endl;
    583 */
    584 //  SetMinimalStringMass2(MinimalStringMass);
    585 }
    586 //*******************************************************************************************************
    587 
    588 void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue)                     // Uzhi
    589 {
    590   MinimalStringMass2=aValue * aValue;
    591 }
    592 //*******************************************************************************************************
    593 
    594 //****************************************************************************************
Note: See TracChangeset for help on using the changeset viewer.