00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 void G4QGSMSplitableHadron::InitParameters()
00052 {
00053
00054 alpha = -0.5;
00055
00056
00057 beta = 2.5;
00058
00059 theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass();
00060
00061
00062
00063 StrangeSuppress = 0.48;
00064 sigmaPt = 0.*GeV;
00065
00066
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
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
00136
00137
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
00144
00145 G4LorentzVector LeftMom(pt, 0.);
00146 G4LorentzVector RightMom;
00147 RightMom.setPx(HadronMom.px() - pt.x());
00148 RightMom.setPy(HadronMom.py() - pt.y());
00149
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
00154 if (Direction) Local2 = -Local2;
00155 G4double RightMinus = 0.5*(Local1 + Local2);
00156 G4double LeftMinus = HadronMom.minus() - RightMinus;
00157
00158
00159 G4double LeftPlus = LeftMom.perp2()/LeftMinus;
00160 G4double RightPlus = HadronMom.plus() - LeftPlus;
00161
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
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
00177 G4double phi, pts;
00178 G4double SumPy = 0.;
00179 G4double SumPx = 0.;
00180 G4ThreeVector Pos = GetPosition();
00181 G4int nSeaPair = GetSoftCollisionCount()-1;
00182
00183
00184
00185
00186
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
00194
00195 G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
00196
00197
00198
00199
00200 G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
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
00221
00222 aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
00223 aParton->SetSpinZ(-firstPartonSpinZ);
00224 aParton->SetColour(-firstPartonColour);
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234 SumPx += aParton->Get4Momentum().px();
00235 SumPy += aParton->Get4Momentum().py();
00236 AntiColor.push_back(aParton);
00237 }
00238
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)
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
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);
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);
00295 SumX += AntiColorX;
00296 if (1. - SumX <= Xmin) break;
00297 }
00298 }
00299 while (1. - SumX <= Xmin);
00300
00301 (*(AntiColor.end()-1))->SetX(1. - SumX);
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
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
00334
00335
00336
00337
00338
00339
00340
00341 Parton2 = new G4Parton(bEnd);
00342 Parton2->SetPosition(GetPosition());
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 Parton2->SetColour(-(Parton1->GetColour()));
00356
00357
00358
00359
00360
00361
00362 if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin())
00363 {
00364 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
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 )
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 }