G4QGSMFragmentation.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // -----------------------------------------------------------------------------
00030 //      GEANT 4 class implementation file
00031 //
00032 //      History: first implementation, Maxim Komogorov, 10-Jul-1998
00033 // -----------------------------------------------------------------------------
00034 #include "G4QGSMFragmentation.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "Randomize.hh"
00037 #include "G4ios.hh"
00038 #include "G4FragmentingString.hh"
00039 #include "G4DiQuarks.hh"
00040 #include "G4Quarks.hh"
00041 
00042 // Class G4QGSMFragmentation 
00043 //****************************************************************************************
00044  
00045 G4QGSMFragmentation::G4QGSMFragmentation() :
00046 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
00047    {
00048    }
00049 
00050 G4QGSMFragmentation::~G4QGSMFragmentation()
00051    {
00052    }
00053 
00054 //----------------------------------------------------------------------------------------------------------
00055 
00056 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
00057 {
00058 //    Can no longer modify Parameters for Fragmentation.
00059         PastInitPhase=true;
00060         
00061 //      check if string has enough mass to fragment...
00062         G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
00063         if ( LeftVector != 0 ) return LeftVector;
00064         
00065         LeftVector = new G4KineticTrackVector;
00066         G4KineticTrackVector * RightVector=new G4KineticTrackVector;
00067 
00068 // this should work but its only a semi deep copy. %GF  G4ExcitedString theStringInCMS(theString);
00069         G4ExcitedString *theStringInCMS=CPExcited(theString);
00070         G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
00071 
00072         G4bool success=false, inner_sucess=true;
00073         G4int attempt=0;
00074         while ( !success && attempt++ < StringLoopInterrupt )
00075         {
00076                 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
00077 
00078                 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
00079                 LeftVector->clear();
00080                 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
00081                 RightVector->clear();
00082                 
00083                 inner_sucess=true;  // set false on failure..
00084                 while (! StopFragmenting(currentString) )
00085                 {  // Split current string into hadron + new string
00086                         G4FragmentingString *newString=0;  // used as output from SplitUp...
00087                         G4KineticTrack * Hadron=Splitup(currentString,newString);
00088                         if ( Hadron != 0 && IsFragmentable(newString)) 
00089                         {
00090                            if ( currentString->GetDecayDirection() > 0 )
00091                                    LeftVector->push_back(Hadron);
00092                            else
00093                                    RightVector->push_back(Hadron);
00094                            delete currentString;
00095                            currentString=newString;
00096                         } else {
00097                          // abandon ... start from the beginning
00098                            if (newString) delete newString;          // Uzhi restore 20.06.08
00099                            if (Hadron)    delete Hadron;
00100                            inner_sucess=false;
00101                            break;
00102                         }
00103                 } 
00104                 // Split current string into 2 final Hadrons
00105                 if ( inner_sucess && 
00106                      SplitLast(currentString,LeftVector, RightVector) ) 
00107                 {
00108                         success=true;
00109                 }
00110                 delete currentString;
00111         }
00112         
00113         delete theStringInCMS;
00114         
00115         if ( ! success )
00116         {
00117                 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
00118                 LeftVector->clear();
00119                 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
00120                 delete RightVector;
00121                 return LeftVector;
00122         }
00123                 
00124         // Join Left- and RightVector into LeftVector in correct order.
00125         while(!RightVector->empty())
00126         {
00127             LeftVector->push_back(RightVector->back());
00128             RightVector->erase(RightVector->end()-1);
00129         }
00130         delete RightVector;
00131 
00132         CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
00133 
00134         G4LorentzRotation toObserverFrame(toCms.inverse());
00135 
00136         for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
00137         {
00138            G4KineticTrack* Hadron = LeftVector->operator[](C1);
00139            G4LorentzVector Momentum = Hadron->Get4Momentum();
00140            Momentum = toObserverFrame*Momentum;
00141            Hadron->Set4Momentum(Momentum);
00142            G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
00143            Momentum = toObserverFrame*Coordinate;
00144            Hadron->SetFormationTime(Momentum.e());
00145            G4ThreeVector aPosition(Momentum.vect());
00146            Hadron->SetPosition(theString.GetPosition()+aPosition);
00147         }
00148         return LeftVector;
00149                 
00150 
00151 
00152 }
00153 
00154 //----------------------------------------------------------------------------------------------------------
00155 
00156 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,  G4ParticleDefinition* , G4double , G4double )
00157 {    
00158   G4double z;    
00159   G4double theA(0), d1, d2, yf;
00160   G4int absCode = std::abs( PartonEncoding );
00161   if (absCode < 10)
00162   { 
00163     if(absCode == 1 || absCode == 2) theA = arho;
00164     else if(absCode == 3) theA = aphi;
00165     else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
00166 
00167     do  
00168     {
00169       z  = zmin + G4UniformRand() * (zmax - zmin);
00170       d1 =  (1. - z);
00171       d2 =  (alft - theA);
00172       yf = std::pow(d1, d2);
00173     } 
00174     while (G4UniformRand() > yf);
00175   }
00176   else
00177   {       
00178     if(absCode == 1103 || absCode == 2101 || 
00179        absCode == 2203 || absCode == 2103)
00180     {
00181       d2 =  (alft - (2.*an - arho));
00182     }
00183     else if(absCode == 3101 || absCode == 3103 ||
00184             absCode == 3201 || absCode == 3203)
00185     {
00186       d2 =  (alft - (2.*ala - arho));
00187     }
00188     else
00189     {
00190       d2 =  (alft - (2.*aksi - arho));
00191     }
00192 
00193     do  
00194     {
00195       z = zmin + G4UniformRand() * (zmax - zmin);
00196       d1 =  (1. - z);
00197       yf = std::pow(d1, d2);
00198     } 
00199     while (G4UniformRand() > yf); 
00200   }
00201   return z;
00202 }
00203 //-----------------------------------------------------------------------------------------
00204 
00205 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
00206                                                   G4FragmentingString * string,    // Uzhi
00207                                                   G4FragmentingString * ) // Uzhi
00208 {
00209        G4double HadronMass = pHadron->GetPDGMass();
00210 
00211        // calculate and assign hadron transverse momentum component HadronPx andHadronPy
00212        G4ThreeVector thePt;
00213        thePt=SampleQuarkPt();
00214        G4ThreeVector HadronPt = thePt +string->DecayPt();
00215        HadronPt.setZ(0);
00216        //...  sample z to define hadron longitudinal momentum and energy
00217        //... but first check the available phase space
00218        G4double DecayQuarkMass2  = sqr(string->GetDecayParton()->GetPDGMass());
00219        G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
00220        if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) ) 
00221           return 0;             // have to start all over!
00222 
00223        //... then compute allowed z region  z_min <= z <= z_max 
00224  
00225        G4double zMin = HadronMass2T/(string->Mass2());
00226        G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
00227        if (zMin >= zMax) return 0;              // have to start all over!
00228         
00229        G4double z = GetLightConeZ(zMin, zMax,
00230                        string->GetDecayParton()->GetPDGEncoding(), pHadron,
00231                        HadronPt.x(), HadronPt.y());      
00232        
00233        //... now compute hadron longitudinal momentum and energy
00234        // longitudinal hadron momentum component HadronPz
00235 
00236         HadronPt.setZ(0.5* string->GetDecayDirection() *
00237                         (z * string->LightConeDecay() - 
00238                          HadronMass2T/(z * string->LightConeDecay())));
00239         G4double HadronE  = 0.5* (z * string->LightConeDecay() + 
00240                                   HadronMass2T/(z * string->LightConeDecay()));
00241 
00242        G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
00243 
00244        return a4Momentum;
00245 }
00246 
00247 
00248 //-----------------------------------------------------------------------------------------
00249 
00250 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
00251                                              G4KineticTrackVector * LeftVector,
00252                                              G4KineticTrackVector * RightVector)
00253 {
00254     //... perform last cluster decay
00255     G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
00256     G4double ResidualMass    =string->Mass(); 
00257     G4double ClusterMassCut = ClusterMass;
00258     G4int cClusterInterrupt = 0;
00259     G4ParticleDefinition * LeftHadron, * RightHadron;
00260     do
00261     {
00262         if (cClusterInterrupt++ >= ClusterLoopInterrupt)
00263         {
00264           return false;
00265         }
00266         G4ParticleDefinition * quark = NULL;
00267         string->SetLeftPartonStable(); // to query quark contents..
00268         if (string->DecayIsQuark() && string->StableIsQuark() ) 
00269         {
00270            //... there are quarks on cluster ends
00271                 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
00272         } else {
00273            //... there is a Diquark on cluster ends
00274                 G4int IsParticle;
00275                 if ( string->StableIsQuark() ) {
00276                   IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 
00277                 } else {
00278                   IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
00279                 }
00280                 pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
00281                 quark = QuarkPair.second;
00282                 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
00283         }
00284         RightHadron = hadronizer->Build(string->GetRightParton(), quark);
00285 
00286        //... repeat procedure, if mass of cluster is too low to produce hadrons
00287        //... ClusterMassCut = 0.15*GeV model parameter
00288         if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
00289         else {ClusterMassCut = ClusterMass;}
00290     } 
00291     while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()  + ClusterMassCut);
00292 
00293     //... compute hadron momenta and energies   
00294     G4LorentzVector  LeftMom, RightMom;
00295     G4ThreeVector    Pos;
00296     Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
00297     LeftMom.boost(ClusterVel);
00298     RightMom.boost(ClusterVel);
00299     LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
00300     RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
00301 
00302     return true;
00303 
00304 }
00305 
00306 //----------------------------------------------------------------------------------------------------------
00307 
00308 G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
00309 {
00310         return sqr(FragmentationMass(string)+MassCut) <
00311                         string->Mass2();
00312 }
00313 
00314 //----------------------------------------------------------------------------------------------------------
00315 
00316 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
00317 {
00318         return
00319          sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
00320          string->Get4Momentum().mag2();
00321 }
00322 
00323 //----------------------------------------------------------------------------------------------------------
00324 
00325 //----------------------------------------------------------------------------------------------------------
00326 
00327 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 
00328     {
00329     G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
00330     G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
00331 
00332     //... sample unit vector       
00333     G4double pz = 1. - 2.*G4UniformRand();  
00334     G4double st     = std::sqrt(1. - pz * pz)*Pabs;
00335     G4double phi    = 2.*pi*G4UniformRand();
00336     G4double px = st*std::cos(phi);
00337     G4double py = st*std::sin(phi);
00338     pz *= Pabs;
00339     
00340     Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
00341     Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
00342 
00343     AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
00344     AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
00345     }
00346 
00347     
00348 //*********************************************************************************************

Generated on Mon May 27 17:49:34 2013 for Geant4 by  doxygen 1.4.7