G4QGSMSplitableHadron.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 #include "G4QGSMSplitableHadron.hh"
00027 #include "G4PhysicalConstants.hh"
00028 #include "G4SystemOfUnits.hh"
00029 #include "G4ParticleTable.hh" 
00030 #include "G4PionPlus.hh"
00031 #include "G4PionMinus.hh"
00032 #include "G4Gamma.hh"
00033 #include "G4PionZero.hh"
00034 #include "G4KaonPlus.hh"
00035 #include "G4KaonMinus.hh"
00036 
00037 // based on prototype by Maxim Komogorov
00038 // Splitting into methods, and centralizing of model parameters HPW Feb 1999
00039 // restructuring HPW Feb 1999
00040 // fixing bug in the sampling of 'x', HPW Feb 1999
00041 // fixing bug in sampling pz, HPW Feb 1999. 
00042 // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
00043 // Using Parton more directly, HPW Feb 1999.
00044 // Shortening the algorithm for sampling x, HPW Feb 1999.
00045 // sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999.
00046 // logic much clearer now. HPW Feb 1999
00047 // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
00048 // Fixing p-t distributions for scattering of nuclei.
00049 // Separating out parameters.
00050 
00051 void G4QGSMSplitableHadron::InitParameters()
00052 {
00053         // changing rapidity distribution for all
00054         alpha = -0.5; // Note that this number is still assumed in the algorithm
00055         // needs to be generalized.
00056         // changing rapidity distribution for projectile like
00057         beta = 2.5;// Note that this number is still assumed in the algorithm
00058         // needs to be generalized.
00059         theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass();
00060         //  theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass();
00061         //  theMinPz = G4PionMinus::PionMinus()->GetPDGMass();
00062         // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
00063         StrangeSuppress = 0.48;
00064         sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7,
00065         // but Maxim's algorithm breaks energy conservation
00066         // to be revised.
00067         widthOfPtSquare = 0.01*GeV*GeV;
00068         Direction = FALSE;
00069         minTransverseMass = 1*keV;
00070 } 
00071 
00072 G4QGSMSplitableHadron::G4QGSMSplitableHadron() 
00073 {
00074         InitParameters();
00075 }
00076 
00077 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection)
00078 :G4VSplitableHadron(aPrimary)
00079 {
00080         InitParameters();
00081         Direction = aDirection;
00082 }
00083 
00084 
00085 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary)
00086 :  G4VSplitableHadron(aPrimary)
00087 {
00088         InitParameters();
00089 }
00090 
00091 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon)
00092 :  G4VSplitableHadron(aNucleon)
00093 {
00094         InitParameters();
00095 }
00096 
00097 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection)
00098 :  G4VSplitableHadron(aNucleon)
00099 {
00100         InitParameters();
00101         Direction = aDirection;
00102 }
00103 
00104 G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){}
00105 
00106 
00107 
00108 //**************************************************************************************************************************
00109 
00110 void G4QGSMSplitableHadron::SplitUp()
00111 {  
00112         if (IsSplit()) return;
00113         Splitting();
00114         if (Color.size()!=0) return;
00115         if (GetSoftCollisionCount() == 0)
00116         {
00117                 DiffractiveSplitUp();
00118         }
00119         else
00120         {
00121                 SoftSplitUp();
00122         }
00123 }
00124 
00125 void G4QGSMSplitableHadron::DiffractiveSplitUp()
00126 {   
00127         // take the particle definitions and get the partons HPW
00128         G4Parton * Left = NULL;
00129         G4Parton * Right = NULL;
00130         GetValenceQuarkFlavors(GetDefinition(), Left, Right);
00131         Left->SetPosition(GetPosition());
00132         Right->SetPosition(GetPosition());
00133 
00134         G4LorentzVector HadronMom = Get4Momentum();
00135         //std::cout << "DSU 1 - "<<HadronMom<<std::endl;
00136 
00137         // momenta of string ends
00138         G4double pt2 = HadronMom.perp2();
00139         G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
00140         G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
00141         G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
00142         if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
00143         //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl;
00144 
00145         G4LorentzVector LeftMom(pt, 0.);
00146         G4LorentzVector RightMom;
00147         RightMom.setPx(HadronMom.px() - pt.x());
00148         RightMom.setPy(HadronMom.py() - pt.y());
00149         //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl;
00150 
00151         G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
00152         G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
00153         //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl;
00154         if (Direction) Local2 = -Local2;
00155         G4double RightMinus   = 0.5*(Local1 + Local2);
00156         G4double LeftMinus = HadronMom.minus() - RightMinus;
00157         //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl;
00158 
00159         G4double LeftPlus  = LeftMom.perp2()/LeftMinus;
00160         G4double RightPlus = HadronMom.plus() - LeftPlus;
00161         //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl;
00162         LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
00163         LeftMom.setE (0.5*(LeftPlus + LeftMinus));
00164         RightMom.setPz(0.5*(RightPlus - RightMinus));
00165         RightMom.setE (0.5*(RightPlus + RightMinus));
00166         //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl;
00167         Left->Set4Momentum(LeftMom);
00168         Right->Set4Momentum(RightMom);
00169         Color.push_back(Left);
00170         AntiColor.push_back(Right);
00171 }
00172 
00173 
00174 void G4QGSMSplitableHadron::SoftSplitUp()
00175 {
00176         //... sample transversal momenta for sea and valence quarks
00177         G4double phi, pts;
00178         G4double SumPy = 0.;
00179         G4double SumPx = 0.;
00180         G4ThreeVector Pos    = GetPosition();
00181         G4int nSeaPair = GetSoftCollisionCount()-1;
00182 
00183         // here the condition,to ensure viability of splitting, also in cases
00184         // where difractive excitation occured together with soft scattering.
00185         //   G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus();
00186         //   G4double Xmin = theMinPz/LightConeMomentum;
00187         G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() );
00188         while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
00189 
00190         G4int aSeaPair;
00191         for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
00192         {
00193                 //  choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
00194 
00195                 G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
00196 
00197                 //  BuildSeaQuark() determines quark spin, isospin and colour
00198                 //  via parton-constructor G4Parton(aPDGCode)
00199 
00200                 G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
00201 
00202                 //              G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
00203 
00204                 //              G4cerr << "Parton 1: "
00205                 //                     << " PDGcode: "  << aPDGCode
00206                 //                     << " - Name: "   << aParton->GetDefinition()->GetParticleName()
00207                 //                     << " - Type: "   << aParton->GetDefinition()->GetParticleType()
00208                 //                     << " - Spin-3: " << aParton->GetSpinZ()
00209                 //                     << " - Colour: " << aParton->GetColour() << G4endl;
00210 
00211                 // save colour a spin-3 for anti-quark
00212 
00213                 G4int firstPartonColour = aParton->GetColour();
00214                 G4double firstPartonSpinZ = aParton->GetSpinZ();
00215 
00216                 SumPx += aParton->Get4Momentum().px();
00217                 SumPy += aParton->Get4Momentum().py();
00218                 Color.push_back(aParton);
00219 
00220                 // create anti-quark
00221 
00222                 aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
00223                 aParton->SetSpinZ(-firstPartonSpinZ);
00224                 aParton->SetColour(-firstPartonColour);
00225 
00226                 //              G4cerr << "Parton 2: "
00227                 //                     << " PDGcode: "  << -aPDGCode
00228                 //                     << " - Name: "   << aParton->GetDefinition()->GetParticleName()
00229                 //                     << " - Type: "   << aParton->GetDefinition()->GetParticleType()
00230                 //                     << " - Spin-3: " << aParton->GetSpinZ()
00231                 //                     << " - Colour: " << aParton->GetColour() << G4endl;
00232                 //              G4cerr << "------------" << G4endl;
00233 
00234                 SumPx += aParton->Get4Momentum().px();
00235                 SumPy += aParton->Get4Momentum().py();
00236                 AntiColor.push_back(aParton);
00237         }
00238         // Valence quark
00239         G4Parton* pColorParton = NULL;
00240         G4Parton* pAntiColorParton = NULL;
00241         GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
00242         G4int ColorEncoding = pColorParton->GetPDGcode();
00243 
00244         pts   =  sigmaPt*std::sqrt(-std::log(G4UniformRand()));
00245         phi   = 2.*pi*G4UniformRand();
00246         G4double Px = pts*std::cos(phi);
00247         G4double Py = pts*std::sin(phi);
00248         SumPx += Px;
00249         SumPy += Py;
00250 
00251         if (ColorEncoding < 0) // use particle definition
00252         {
00253                 G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0);
00254                 pColorParton->Set4Momentum(ColorMom);
00255                 G4LorentzVector AntiColorMom(Px, Py, 0, 0);
00256                 pAntiColorParton->Set4Momentum(AntiColorMom);
00257         }
00258         else
00259         {
00260                 G4LorentzVector ColorMom(Px, Py, 0, 0);
00261                 pColorParton->Set4Momentum(ColorMom);
00262                 G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0);
00263                 pAntiColorParton->Set4Momentum(AntiColorMom);
00264         }
00265         Color.push_back(pColorParton);
00266         AntiColor.push_back(pAntiColorParton);
00267 
00268         // Sample X
00269         G4int nAttempt = 0;
00270         G4double SumX = 0;
00271         G4double aBeta = beta;
00272         G4double ColorX, AntiColorX;
00273         if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.;
00274         if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.;
00275         if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.;
00276         if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.;
00277         if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;
00278         if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;
00279         do
00280         {
00281                 SumX = 0;
00282                 nAttempt++;
00283                 G4int    NumberOfUnsampledSeaQuarks = 2*nSeaPair;
00284                 ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
00285                 Color.back()->SetX(SumX = ColorX);// this is the valenz quark.
00286                 for(G4int aPair = 0; aPair < nSeaPair; aPair++)
00287                 {
00288                         NumberOfUnsampledSeaQuarks--;
00289                         ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
00290                         Color[aPair]->SetX(ColorX);
00291                         SumX += ColorX;
00292                         NumberOfUnsampledSeaQuarks--;
00293                         AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
00294                         AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons
00295                         SumX += AntiColorX;
00296                         if (1. - SumX <= Xmin)  break;
00297                 }
00298         }
00299         while (1. - SumX <= Xmin);
00300 
00301         (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum
00302         G4double lightCone  = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus());
00303         G4double lightCone2 = ((!Direction) ? Get4Momentum().plus() : Get4Momentum().minus());
00304         for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++)
00305         {
00306                 G4Parton* aParton = Color[aSeaPair];
00307                 aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
00308 
00309                 aParton = AntiColor[aSeaPair];
00310                 aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
00311         }
00312         return;
00313 } 
00314 
00315 void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2)
00316 {
00317         // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
00318         G4int aEnd;
00319         G4int bEnd;
00320         G4int HadronEncoding = aPart->GetPDGEncoding();
00321         if (aPart->GetBaryonNumber() == 0)
00322         {
00323                 theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
00324         }
00325         else
00326         {
00327                 theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
00328         }
00329 
00330         Parton1 = new G4Parton(aEnd);
00331         Parton1->SetPosition(GetPosition());
00332 
00333         //      G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
00334         //      G4cerr << "Parton 1: "
00335         //             << " PDGcode: "  << aEnd
00336         //             << " - Name: "   << Parton1->GetDefinition()->GetParticleName()
00337         //             << " - Type: "   << Parton1->GetDefinition()->GetParticleType()
00338         //             << " - Spin-3: " << Parton1->GetSpinZ()
00339         //             << " - Colour: " << Parton1->GetColour() << G4endl;
00340 
00341         Parton2 = new G4Parton(bEnd);
00342         Parton2->SetPosition(GetPosition());
00343 
00344         //      G4cerr << "Parton 2: "
00345         //             << " PDGcode: "  << bEnd
00346         //             << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
00347         //             << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
00348         //             << " - Spin-3: " << Parton2->GetSpinZ()
00349         //             << " - Colour: " << Parton2->GetColour() << G4endl;
00350         //      G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
00351 
00352         // colour of parton 1 choosen at random by G4Parton(aEnd)
00353         // colour of parton 2 is the opposite:
00354 
00355         Parton2->SetColour(-(Parton1->GetColour()));
00356 
00357         // isospin-3 of both partons is handled by G4Parton(PDGCode)
00358 
00359         // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
00360         // spin-3 of parton 2 may be constrained by spin of original particle:
00361 
00362         if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin())
00363         {
00364                 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));    
00365         }
00366 
00367         //      G4cerr << "Parton 2: "
00368         //             << " PDGcode: "  << bEnd
00369         //             << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
00370         //             << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
00371         //             << " - Spin-3: " << Parton2->GetSpinZ()
00372         //             << " - Colour: " << Parton2->GetColour() << G4endl;
00373         //      G4cerr << "------------" << G4endl;
00374 
00375 }
00376 
00377 
00378 G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
00379 {
00380         G4double R;
00381         while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) {;}
00382         R = std::sqrt(R);
00383         G4double phi = twopi*G4UniformRand();
00384         return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);
00385 }
00386 
00387 G4Parton * G4QGSMSplitableHadron::
00388 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
00389 {
00390         if (isAntiQuark) aPDGCode*=-1;
00391         G4Parton* result = new G4Parton(aPDGCode);
00392         result->SetPosition(GetPosition());
00393         G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
00394         G4LorentzVector a4Momentum(aPtVector, 0);
00395         result->Set4Momentum(a4Momentum);
00396         return result;
00397 }
00398 
00399 G4double G4QGSMSplitableHadron::
00400 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
00401 {
00402         G4double result;
00403         G4double x1, x2;
00404         G4double ymax = 0;
00405         for(G4int ii=1; ii<100; ii++)
00406         {
00407                 G4double y = std::pow(1./G4double(ii), alpha);
00408                 y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
00409                 y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
00410                 if(y>ymax) ymax = y;
00411         }
00412         G4double y;
00413         G4double xMax=1-(totalSea+1)*anXmin;
00414         if(anXmin > xMax)
00415         {
00416                 G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl;
00417                 throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
00418         }
00419         do
00420         {
00421                 x1 = CLHEP::RandFlat::shoot(anXmin, xMax);
00422                 y = std::pow(x1, alpha);
00423                 y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
00424                 y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
00425                 x2 = ymax*G4UniformRand();
00426         }
00427         while(x2>y);
00428         result = x1;
00429         return result;
00430 }

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