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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "globals.hh"
00047 #include "Randomize.hh"
00048 #include "G4PhysicalConstants.hh"
00049 #include "G4SystemOfUnits.hh"
00050
00051 #include "G4DiffractiveExcitation.hh"
00052 #include "G4FTFParameters.hh"
00053 #include "G4ElasticHNScattering.hh"
00054
00055 #include "G4LorentzRotation.hh"
00056 #include "G4RotationMatrix.hh"
00057 #include "G4ThreeVector.hh"
00058 #include "G4ParticleDefinition.hh"
00059 #include "G4VSplitableHadron.hh"
00060 #include "G4ExcitedString.hh"
00061 #include "G4ParticleTable.hh"
00062 #include "G4Neutron.hh"
00063 #include "G4ParticleDefinition.hh"
00064
00065
00066
00067
00068 G4DiffractiveExcitation::G4DiffractiveExcitation()
00069 {
00070 }
00071
00072
00073 G4bool G4DiffractiveExcitation::
00074 ExciteParticipants(G4VSplitableHadron *projectile,
00075 G4VSplitableHadron *target,
00076 G4FTFParameters *theParameters,
00077 G4ElasticHNScattering *theElastic) const
00078 {
00079
00080
00081 G4LorentzVector Pprojectile=projectile->Get4Momentum();
00082
00083 if(Pprojectile.z() < 0.)
00084 {
00085 target->SetStatus(2);
00086 return false;
00087 }
00088
00089 G4double ProjectileRapidity = Pprojectile.rapidity();
00090
00091 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
00092 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode);
00093
00094 G4bool PutOnMassShell(false);
00095
00096 G4double M0projectile = Pprojectile.mag();
00097
00098
00099 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
00100 {
00101 PutOnMassShell=true;
00102 M0projectile=projectile->GetDefinition()->GetPDGMass();
00103 }
00104
00105 G4double M0projectile2 = M0projectile * M0projectile;
00106
00107 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
00108 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
00109 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
00110
00111
00112 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
00113 G4int absTargetPDGcode=std::abs(TargetPDGcode);
00114
00115
00116 G4LorentzVector Ptarget=target->Get4Momentum();
00117
00118
00119 G4double M0target = Ptarget.mag();
00120
00121
00122
00123 if(M0target < target->GetDefinition()->GetPDGMass())
00124 {
00125 PutOnMassShell=true;
00126 M0target=target->GetDefinition()->GetPDGMass();
00127 }
00128
00129 G4double M0target2 = M0target * M0target;
00130
00131 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();
00132 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();
00133 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
00134
00135 G4double AveragePt2=theParameters->GetAveragePt2();
00136 G4double ProbLogDistr=theParameters->GetProbLogDistr();
00137
00138
00139
00140
00141
00142 G4double SumMasses=M0projectile+M0target+220.*MeV;
00143
00144
00145 G4LorentzVector Psum;
00146 Psum=Pprojectile+Ptarget;
00147 G4double S=Psum.mag2();
00148
00149
00150
00151
00152 G4LorentzRotation toCms(-1*Psum.boostVector());
00153
00154 G4LorentzVector Ptmp=toCms*Pprojectile;
00155
00156 if ( Ptmp.pz() <= 0. )
00157 {
00158 target->SetStatus(2);
00159
00160 return false;
00161 }
00162
00163 toCms.rotateZ(-1*Ptmp.phi());
00164 toCms.rotateY(-1*Ptmp.theta());
00165
00166 G4LorentzRotation toLab(toCms.inverse());
00167
00168 Pprojectile.transform(toCms);
00169 Ptarget.transform(toCms);
00170
00171 G4double PZcms2, PZcms;
00172
00173 G4double SqrtS=std::sqrt(S);
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 if(absProjectilePDGcode > 1000 && (SqrtS < SumMasses))
00184 {target->SetStatus(2); return false;}
00185
00186
00187
00188
00189
00190
00191
00192 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) &&
00193 (SqrtS < SumMasses))
00194 {target->SetStatus(2); return false;}
00195
00196
00197
00198
00199 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) ||
00200 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
00201 (absProjectilePDGcode == 310)) && (SqrtS < SumMasses))
00202
00203 {target->SetStatus(2); return false;}
00204
00205
00206
00207 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
00208 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
00209 /4./S;
00210
00211 if(PZcms2 < 0)
00212 {target->SetStatus(2); return false;}
00213
00214
00215 PZcms = std::sqrt(PZcms2);
00216
00217 if(PutOnMassShell)
00218 {
00219 if(Pprojectile.z() > 0.)
00220 {
00221 Pprojectile.setPz( PZcms);
00222 Ptarget.setPz( -PZcms);
00223 }
00224 else
00225 {
00226 Pprojectile.setPz(-PZcms);
00227 Ptarget.setPz( PZcms);
00228 };
00229
00230 Pprojectile.setE(std::sqrt(M0projectile2 +
00231 Pprojectile.x()*Pprojectile.x()+
00232 Pprojectile.y()*Pprojectile.y()+
00233 PZcms2));
00234 Ptarget.setE(std::sqrt(M0target2 +
00235 Ptarget.x()*Ptarget.x()+
00236 Ptarget.y()*Ptarget.y()+
00237 PZcms2));
00238 }
00239
00240
00241
00242
00243 G4double maxPtSquare;
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange();
00262 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange();
00263
00264 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
00265
00266
00267
00268 G4double DeltaMass=
00269 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 if(G4UniformRand() < MagQuarkExchange*
00281 std::exp(-SlopeQuarkExchange*ProjectileRapidity))
00282 {
00283
00284
00285
00286 G4int NewProjCode(0), NewTargCode(0);
00287
00288 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
00289
00290
00291 if(absProjectilePDGcode < 1000 )
00292 {
00293 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2);
00294 } else
00295 {
00296 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
00297 }
00298
00299
00300 G4int TargQ1(0), TargQ2(0), TargQ3(0);
00301 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3);
00302
00303
00304
00305
00306 G4int ProjExchangeQ(0);
00307 G4int TargExchangeQ(0);
00308
00309 if(absProjectilePDGcode < 1000 )
00310 {
00311
00312 if(ProjQ1 > 0 )
00313 {
00314 G4int Navailable=0;
00315 ProjExchangeQ = ProjQ1;
00316 if(ProjExchangeQ != TargQ1) Navailable++;
00317 if(ProjExchangeQ != TargQ2) Navailable++;
00318 if(ProjExchangeQ != TargQ3) Navailable++;
00319
00320 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
00321
00322 Navailable=0;
00323 if(ProjExchangeQ != TargQ1)
00324 {
00325 Navailable++;
00326 if(Navailable == Nsampled)
00327 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00328 }
00329
00330 if(ProjExchangeQ != TargQ2)
00331 {
00332 Navailable++;
00333 if(Navailable == Nsampled)
00334 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00335 }
00336
00337 if(ProjExchangeQ != TargQ3)
00338 {
00339 Navailable++;
00340 if(Navailable == Nsampled)
00341 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00342 }
00343 } else
00344 {
00345 G4int Navailable=0;
00346 ProjExchangeQ = ProjQ2;
00347 if(ProjExchangeQ != TargQ1) Navailable++;
00348 if(ProjExchangeQ != TargQ2) Navailable++;
00349 if(ProjExchangeQ != TargQ3) Navailable++;
00350
00351 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
00352
00353 Navailable=0;
00354 if(ProjExchangeQ != TargQ1)
00355 {
00356 Navailable++;
00357 if(Navailable == Nsampled)
00358 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00359 }
00360
00361 if(ProjExchangeQ != TargQ2)
00362 {
00363 Navailable++;
00364 if(Navailable == Nsampled)
00365 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00366 }
00367
00368 if(ProjExchangeQ != TargQ3)
00369 {
00370 Navailable++;
00371 if(Navailable == Nsampled)
00372 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00373 }
00374 }
00375
00376
00377
00378
00379
00380 G4int aProjQ1=std::abs(ProjQ1);
00381 G4int aProjQ2=std::abs(ProjQ2);
00382 if(aProjQ1 == aProjQ2) {NewProjCode = 111;}
00383 else
00384 {
00385 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
00386 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
00387
00388 }
00389
00390 G4bool ProjExcited=false;
00391
00392 if(G4UniformRand() < 0.5)
00393 {
00394 NewProjCode +=2;
00395 ProjExcited=true;
00396 }
00397 if(aProjQ1 != aProjQ2) NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode);
00398
00399
00400 G4ParticleDefinition* TestParticle=0;
00401 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
00402
00403
00404 if(TestParticle)
00405 {
00406 G4double MtestPart=
00407 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
00408
00409
00410
00411
00412
00413
00414
00415 if(MtestPart > M0projectile)
00416 {M0projectile = MtestPart;}
00417 else
00418 {
00419 if(std::abs(M0projectile - projectile->GetDefinition()->GetPDGMass()) < 140.*MeV)
00420 {M0projectile = MtestPart;}
00421 }
00422
00423 M0projectile2 = M0projectile * M0projectile;
00424
00425 ProjectileDiffStateMinMass =M0projectile+210.*MeV;
00426 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV;
00427 } else
00428 {return false;}
00429
00430
00431 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
00432
00433
00434
00435
00436
00437
00438 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
00439 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;
00440 ProjExcited=true;}
00441 else if( target->GetDefinition()->GetPDGiIsospin() == 3 )
00442 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;
00443 ProjExcited=true;}
00444 else {}
00445 }
00446
00447
00448 else if((!ProjExcited) &&
00449 (G4UniformRand() < DeltaProbAtQuarkExchange) &&
00450 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;}
00451
00452 else {}
00453
00454
00455
00456
00457 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
00458
00459 if(TestParticle)
00460 {
00461 G4double MtestPart=
00462 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
00463
00464 if(MtestPart > M0target)
00465 {M0target=MtestPart;}
00466 else
00467 {
00468 if(std::abs(M0target - target->GetDefinition()->GetPDGMass()) < 140.*MeV)
00469 {M0target=MtestPart;}
00470 }
00471
00472 TargetDiffStateMinMass =M0target+220.*MeV;
00473 TargetNonDiffStateMinMass=M0target+220.*MeV;
00474 } else
00475 {return false;}
00476 } else
00477 {
00478
00479 G4double Same=theParameters->GetProbOfSameQuarkExchange();
00480 G4bool ProjDeltaHasCreated(false);
00481 G4bool TargDeltaHasCreated(false);
00482
00483 G4double Ksi=G4UniformRand();
00484 if(G4UniformRand() < 0.5)
00485 {
00486 if( Ksi < 0.333333 )
00487 {ProjExchangeQ = ProjQ1;}
00488 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00489 {ProjExchangeQ = ProjQ2;}
00490 else
00491 {ProjExchangeQ = ProjQ3;}
00492
00493
00494 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same))
00495 {
00496 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00497 } else
00498 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same))
00499 {
00500 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00501 } else
00502 {
00503 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00504 }
00505
00506
00507
00508 if( Ksi < 0.333333 )
00509 {ProjQ1=ProjExchangeQ;}
00510 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00511 {ProjQ2=ProjExchangeQ;}
00512 else
00513 {ProjQ3=ProjExchangeQ;}
00514
00515 } else
00516 {
00517 if( Ksi < 0.333333 )
00518 {TargExchangeQ = TargQ1;}
00519 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00520 {TargExchangeQ = TargQ2;}
00521 else
00522 {TargExchangeQ = TargQ3;}
00523
00524 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same))
00525 {
00526 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00527 } else
00528 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same))
00529 {
00530 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00531 } else
00532 {
00533 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00534 }
00535
00536 if( Ksi < 0.333333 )
00537 {TargQ1=TargExchangeQ;}
00538 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00539 {TargQ2=TargExchangeQ;}
00540 else
00541 {TargQ3=TargExchangeQ;}
00542
00543 }
00544
00545 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3);
00546
00547 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
00548 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)
00549 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
00550 {NewProjCode +=2; ProjDeltaHasCreated=true;}
00551 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
00552 }
00553 else
00554 {
00555 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
00556 {NewProjCode +=2; ProjDeltaHasCreated=true;}
00557 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
00558 }
00559
00560
00561
00562
00563
00564
00565
00566
00567 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}
00578 else if(target->GetDefinition()->GetPDGiIsospin() == 3)
00579 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
00580 {NewTargCode +=2; TargDeltaHasCreated=true;}
00581 else {NewTargCode +=0; TargDeltaHasCreated=false;}
00582 }
00583 else
00584 {
00585 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
00586 {NewTargCode +=2; TargDeltaHasCreated=true;}
00587 else {NewTargCode +=0; TargDeltaHasCreated=false;}
00588 }
00589
00590
00591
00592
00593
00594
00595
00596
00597 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
00598 {
00599 }
00600
00601 if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;}
00602 if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;}
00603 if(ProjDeltaHasCreated)
00604 {
00605 G4double MtestPart=
00606 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
00607
00608 if(MtestPart >= M0projectile)
00609 {
00610 M0projectile = MtestPart;
00611 M0projectile2 = M0projectile * M0projectile;
00612 }
00613
00614 ProjectileDiffStateMinMass =M0projectile+210.*MeV;
00615 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV;
00616 }
00617
00618
00619
00620 if(TargDeltaHasCreated)
00621 {
00622 G4double MtestPart=
00623 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
00624
00625 if(MtestPart >=M0target)
00626 {
00627 M0target=MtestPart;
00628 M0target2 = M0target * M0target;
00629 }
00630
00631 TargetDiffStateMinMass =M0target+210.*MeV;
00632 TargetNonDiffStateMinMass=M0target+210.*MeV;
00633 }
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
00644 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
00645 /4./S;
00646
00647 if(PZcms2 < 0) {return false;}
00648
00649 projectile->SetDefinition(
00650 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
00651
00652 target->SetDefinition(
00653 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));
00654
00655 PZcms = std::sqrt(PZcms2);
00656
00657 Pprojectile.setPz( PZcms);
00658 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
00659
00660 Ptarget.setPz( -PZcms);
00661 Ptarget.setE(std::sqrt(M0target2+PZcms2));
00662
00663
00664
00665 if(absProjectilePDGcode < 1000)
00666 {
00667 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
00668 Wexcit=0.;
00669 if(G4UniformRand() > Wexcit)
00670 {
00671
00672 Pprojectile.transform(toLab);
00673 Ptarget.transform(toLab);
00674
00675 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00676 projectile->SetPosition(target->GetPosition());
00677
00678 projectile->Set4Momentum(Pprojectile);
00679 target->Set4Momentum(Ptarget);
00680
00681 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
00682
00683 return Result;
00684 }
00685 } else
00686 {
00687 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
00688
00689 if(G4UniformRand() > Wexcit)
00690 {
00691
00692 Pprojectile.transform(toLab);
00693 Ptarget.transform(toLab);
00694
00695 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00696 projectile->SetPosition(target->GetPosition());
00697
00698 projectile->Set4Momentum(Pprojectile);
00699 target->Set4Momentum(Ptarget);
00700
00701 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
00702 return Result;
00703 }
00704 }
00705
00706 }
00707
00708
00709 G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction;
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732 if(ProbOfDiffraction!=0.)
00733 {
00734 ProbProjectileDiffraction/=ProbOfDiffraction;
00735 }
00736 else
00737 {
00738 ProbProjectileDiffraction=0.;
00739 }
00740
00741
00742
00743 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass *
00744 ProjectileDiffStateMinMass;
00745 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
00746 ProjectileNonDiffStateMinMass;
00747
00748 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass *
00749 TargetDiffStateMinMass;
00750 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass *
00751 TargetNonDiffStateMinMass;
00752
00753 G4double Pt2;
00754 G4double ProjMassT2, ProjMassT;
00755 G4double TargMassT2, TargMassT;
00756 G4double PMinusMin, PMinusMax;
00757
00758 G4double TPlusMin , TPlusMax;
00759 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
00760
00761 G4LorentzVector Qmomentum;
00762 G4double Qminus, Qplus;
00763
00764 G4int whilecount=0;
00765
00766
00767
00768
00769 if(G4UniformRand() < ProbOfDiffraction)
00770 {
00771 if(G4UniformRand() < ProbProjectileDiffraction)
00772 {
00773
00774
00775 do {
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788 if (whilecount > 1000 )
00789 {
00790 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00791 target->SetStatus(2); return false;
00792 };
00793
00794
00795 ProjMassT2=ProjectileDiffStateMinMass2;
00796 ProjMassT =ProjectileDiffStateMinMass;
00797
00798 TargMassT2=M0target2;
00799 TargMassT =M0target;
00800
00801 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00802 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00803 /4./S;
00804
00805
00806 if(PZcms2 < 0 )
00807 {
00808 target->SetStatus(2);
00809 return false;
00810 }
00811 maxPtSquare=PZcms2;
00812
00813 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00814 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00815
00816 ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
00817 ProjMassT =std::sqrt(ProjMassT2);
00818
00819 TargMassT2=M0target2+Pt2;
00820 TargMassT =std::sqrt(TargMassT2);
00821
00822
00823
00824 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00825 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00826 /4./S;
00827
00828 if(PZcms2 < 0 ) continue;
00829 PZcms =std::sqrt(PZcms2);
00830
00831 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
00832 PMinusMax=SqrtS-TargMassT;
00833
00834 PMinusNew=ChooseP(PMinusMin, PMinusMax);
00835
00836
00837
00838 TMinusNew=SqrtS-PMinusNew;
00839 Qminus=Ptarget.minus()-TMinusNew;
00840 TPlusNew=TargMassT2/TMinusNew;
00841 Qplus=Ptarget.plus()-TPlusNew;
00842
00843 Qmomentum.setPz( (Qplus-Qminus)/2 );
00844 Qmomentum.setE( (Qplus+Qminus)/2 );
00845
00846 } while ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2);
00847
00848
00849 projectile->SetStatus(1*projectile->GetStatus());
00850 }
00851 else
00852 {
00853
00854
00855 do {
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868 if (whilecount > 1000 )
00869 {
00870 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00871 target->SetStatus(2); return false;
00872 };
00873
00874
00875 ProjMassT2=M0projectile2;
00876 ProjMassT =M0projectile;
00877
00878 TargMassT2=TargetDiffStateMinMass2;
00879 TargMassT =TargetDiffStateMinMass;
00880
00881 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00882 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00883 /4./S;
00884
00885
00886 if(PZcms2 < 0 )
00887 {
00888 target->SetStatus(2);
00889 return false;
00890 }
00891 maxPtSquare=PZcms2;
00892
00893 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00894
00895
00896 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00897
00898 ProjMassT2=M0projectile2+Pt2;
00899 ProjMassT =std::sqrt(ProjMassT2);
00900
00901 TargMassT2=TargetDiffStateMinMass2+Pt2;
00902 TargMassT =std::sqrt(TargMassT2);
00903
00904 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00905 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00906 /4./S;
00907
00908
00909 if(PZcms2 < 0 ) continue;
00910 PZcms =std::sqrt(PZcms2);
00911
00912 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
00913 TPlusMax=SqrtS-ProjMassT;
00914
00915 TPlusNew=ChooseP(TPlusMin, TPlusMax);
00916
00917
00918
00919
00920 PPlusNew=SqrtS-TPlusNew;
00921 Qplus=PPlusNew-Pprojectile.plus();
00922 PMinusNew=ProjMassT2/PPlusNew;
00923 Qminus=PMinusNew-Pprojectile.minus();
00924
00925 Qmomentum.setPz( (Qplus-Qminus)/2 );
00926 Qmomentum.setE( (Qplus+Qminus)/2 );
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938 } while ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2);
00939
00940
00941
00942
00943 target->SetStatus(1*target->GetStatus());
00944 }
00945 }
00946 else
00947 {
00948
00949
00950 do {
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963 if (whilecount > 1000 )
00964 {
00965 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00966 target->SetStatus(2); return false;
00967 };
00968
00969 ProjMassT2=ProjectileNonDiffStateMinMass2;
00970 ProjMassT =ProjectileNonDiffStateMinMass;
00971
00972 TargMassT2=TargetNonDiffStateMinMass2;
00973 TargMassT =TargetNonDiffStateMinMass;
00974
00975 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00976 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00977 /4./S;
00978
00979 if(PZcms2 < 0 )
00980 {
00981 target->SetStatus(2);
00982 return false;
00983 }
00984 maxPtSquare=PZcms2;
00985 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00986 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00987
00988 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
00989 ProjMassT =std::sqrt(ProjMassT2);
00990
00991 TargMassT2=TargetNonDiffStateMinMass2+Pt2;
00992 TargMassT =std::sqrt(TargMassT2);
00993
00994 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00995 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00996 /4./S;
00997
00998
00999 if(PZcms2 < 0 ) continue;
01000 PZcms =std::sqrt(PZcms2);
01001
01002 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
01003 PMinusMax=SqrtS-TargMassT;
01004
01005 if(G4UniformRand() < ProbLogDistr)
01006 { PMinusNew=ChooseP(PMinusMin, PMinusMax);}
01007 else {PMinusNew=(PMinusMax-PMinusMin)*G4UniformRand() + PMinusMin;}
01008 Qminus=PMinusNew-Pprojectile.minus();
01009
01010 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
01011 TPlusMax=SqrtS-PMinusNew;
01012
01013
01014 if(G4UniformRand() < 0.5)
01015 { TPlusNew=ChooseP(TPlusMin, TPlusMax);}
01016 else {TPlusNew=(TPlusMax-TPlusMin)*G4UniformRand() +TPlusMin;}
01017
01018 Qplus=-(TPlusNew-Ptarget.plus());
01019
01020 Qmomentum.setPz( (Qplus-Qminus)/2 );
01021 Qmomentum.setE( (Qplus+Qminus)/2 );
01022
01023
01024
01025
01026
01027 } while (
01028 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) ||
01029 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 ));
01030
01031 projectile->SetStatus(0*projectile->GetStatus());
01032 target->SetStatus(0*target->GetStatus());
01033 }
01034
01035 Pprojectile += Qmomentum;
01036 Ptarget -= Qmomentum;
01037
01038
01039
01040
01041 Pprojectile.transform(toLab);
01042 Ptarget.transform(toLab);
01043
01044
01045 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
01046 projectile->SetPosition(target->GetPosition());
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066 projectile->Set4Momentum(Pprojectile);
01067 target->Set4Momentum(Ptarget);
01068
01069 projectile->IncrementCollisionCount(1);
01070 target->IncrementCollisionCount(1);
01071
01072 return true;
01073 }
01074
01075
01076 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron,
01077 G4bool isProjectile,
01078 G4ExcitedString * &FirstString,
01079 G4ExcitedString * &SecondString,
01080 G4FTFParameters *theParameters) const
01081 {
01082
01083
01084
01085
01086
01087 hadron->SplitUp();
01088 G4Parton *start= hadron->GetNextParton();
01089
01090 if ( start==NULL)
01091 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
01092 FirstString=0; SecondString=0;
01093 return;
01094 }
01095 G4Parton *end = hadron->GetNextParton();
01096 if ( end==NULL)
01097 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
01098 FirstString=0; SecondString=0;
01099 return;
01100 }
01101
01102
01103 G4LorentzVector Phadron=hadron->Get4Momentum();
01104
01105 G4LorentzVector Pstart(0.,0.,0.,0.);
01106 G4LorentzVector Pend(0.,0.,0.,0.);
01107 G4LorentzVector Pkink(0.,0.,0.,0.);
01108 G4LorentzVector PkinkQ1(0.,0.,0.,0.);
01109 G4LorentzVector PkinkQ2(0.,0.,0.,0.);
01110
01111 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
01112 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding());
01113
01114
01115
01116 G4double Wmin(0.);
01117 if(isProjectile)
01118 {
01119 Wmin=theParameters->GetProjMinDiffMass();
01120 } else
01121 {
01122 Wmin=theParameters->GetTarMinDiffMass();
01123 }
01124
01125 G4double W = hadron->Get4Momentum().mag();
01126
01127 G4double W2=W*W;
01128
01129 G4double Pt(0.), x1(0.), x3(0.);
01130 G4bool Kink=false;
01131
01132 if(!(((start->GetDefinition()->GetParticleSubType() == "di_quark") &&
01133 ( end->GetDefinition()->GetParticleSubType() == "di_quark") ) ||
01134 ((start->GetDefinition()->GetParticleSubType() == "quark") &&
01135 ( end->GetDefinition()->GetParticleSubType() == "quark") )))
01136 {
01137
01138
01139
01140 if(W > Wmin)
01141 {
01142 if(hadron->GetStatus() == 0)
01143 {
01144 G4double Pt2kink=theParameters->GetPt2Kink();
01145 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
01146 }
01147 else
01148 {Pt=0.;}
01149
01150 if(Pt > 500.*MeV)
01151 {
01152 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
01153 G4double Y=Ymax*(1.- 2.*G4UniformRand());
01154
01155 x1=1.-Pt/W*std::exp( Y);
01156 x3=1.-Pt/W*std::exp(-Y);
01157
01158
01159 G4double Mass_startQ = 650.*MeV;
01160 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV;
01161 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV;
01162 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
01163
01164 G4double Mass_endQ = 650.*MeV;
01165 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV;
01166 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV;
01167 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
01168
01169 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ;
01170 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
01171
01172 G4double P2_2=sqr((2.-x1-x3)*W/2.);
01173
01174 if((P2_1 <= 0.) || (P2_3 <= 0.))
01175 { Kink=false;}
01176 else
01177 {
01178 G4double P_1=std::sqrt(P2_1);
01179 G4double P_2=std::sqrt(P2_2);
01180 G4double P_3=std::sqrt(P2_3);
01181
01182 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
01183 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
01184
01185
01186 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.))
01187 { Kink=false;}
01188 else
01189 {
01190 Kink=true;
01191 Pt=P_2*std::sqrt(1.-CosT12*CosT12);
01192 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13);
01193 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1);
01194 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12);
01195 Pstart.setE(x3*W/2.);
01196 Pkink.setE(Pkink.vect().mag());
01197 Pend.setE(x1*W/2.);
01198
01199 G4double XkQ=GetQuarkFractionOfKink(0.,1.);
01200 if(Pkink.getZ() > 0.)
01201 {
01202 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
01203 } else {
01204 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
01205 }
01206
01207 PkinkQ2=Pkink - PkinkQ1;
01208
01209
01210 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
01211 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
01212 G4double Psi=std::acos(Cos2Psi);
01213
01214 G4LorentzRotation Rotate;
01215 if(isProjectile) {Rotate.rotateY(Psi);}
01216 else {Rotate.rotateY(pi-Psi);}
01217 Rotate.rotateZ(twopi*G4UniformRand());
01218
01219 Pstart*=Rotate;
01220 Pkink*=Rotate;
01221 PkinkQ1*=Rotate;
01222 PkinkQ2*=Rotate;
01223 Pend*=Rotate;
01224
01225 }
01226 }
01227 }
01228 }
01229 }
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242 if(Kink)
01243 {
01244
01245 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
01246 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
01247
01248 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
01249 for(unsigned int Iq=0; Iq <3; Iq++)
01250 {
01251
01252
01253 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
01254
01255
01256 G4Parton * Gquark = new G4Parton(QuarkInGluon);
01257 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
01258
01259
01260
01261 G4LorentzRotation toCMS(-1*Phadron.boostVector());
01262
01263 G4LorentzRotation toLab(toCMS.inverse());
01264
01265
01266 Pstart.transform(toLab); start->Set4Momentum(Pstart);
01267 PkinkQ1.transform(toLab);
01268 PkinkQ2.transform(toLab);
01269 Pend.transform(toLab); end->Set4Momentum(Pend);
01270
01271
01272
01273
01274
01275
01276 G4int absPDGcode=1500;
01277 if((start->GetDefinition()->GetParticleSubType() == "quark") &&
01278 ( end->GetDefinition()->GetParticleSubType() == "quark") )
01279 absPDGcode=110;
01280
01281
01282
01283 if(absPDGcode < 1000)
01284 {
01285 if ( isProjectile )
01286 {
01287 if(end->GetDefinition()->GetPDGEncoding() > 0 )
01288 {
01289 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
01290 SecondString= new G4ExcitedString(Gquark,start ,+1);
01291 Ganti_quark->Set4Momentum(PkinkQ1);
01292 Gquark->Set4Momentum(PkinkQ2);
01293
01294 } else
01295 {
01296 FirstString = new G4ExcitedString(end ,Gquark, +1);
01297 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
01298 Gquark->Set4Momentum(PkinkQ1);
01299 Ganti_quark->Set4Momentum(PkinkQ2);
01300
01301 }
01302 } else {
01303 if(end->GetDefinition()->GetPDGEncoding() > 0 )
01304 {
01305 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
01306 SecondString= new G4ExcitedString(start ,Gquark,-1);
01307 Ganti_quark->Set4Momentum(PkinkQ2);
01308 Gquark->Set4Momentum(PkinkQ1);
01309
01310 } else
01311 {
01312 FirstString = new G4ExcitedString(Gquark,end ,-1);
01313 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
01314 Gquark->Set4Momentum(PkinkQ2);
01315 Ganti_quark->Set4Momentum(PkinkQ1);
01316
01317 }
01318 }
01319 } else
01320 {
01321 if ( isProjectile )
01322 {
01323 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
01324 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
01325 {
01326 FirstString = new G4ExcitedString(end ,Gquark, +1);
01327 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
01328 Gquark->Set4Momentum(PkinkQ1);
01329 Ganti_quark->Set4Momentum(PkinkQ2);
01330
01331 } else
01332 {
01333 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
01334 SecondString= new G4ExcitedString(Gquark,start ,+1);
01335 Ganti_quark->Set4Momentum(PkinkQ1);
01336 Gquark->Set4Momentum(PkinkQ2);
01337
01338 }
01339 } else {
01340
01341 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
01342 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
01343 {
01344 FirstString = new G4ExcitedString(Gquark,end ,-1);
01345
01346 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
01347 Gquark->Set4Momentum(PkinkQ1);
01348 Ganti_quark->Set4Momentum(PkinkQ2);
01349
01350 } else
01351 {
01352 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
01353 SecondString= new G4ExcitedString(start ,Gquark,-1);
01354 Gquark->Set4Momentum(PkinkQ2);
01355 Ganti_quark->Set4Momentum(PkinkQ1);
01356
01357 }
01358 }
01359 }
01360
01361 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01362 FirstString->SetPosition(hadron->GetPosition());
01363
01364 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01365 SecondString->SetPosition(hadron->GetPosition());
01366
01367
01368 } else
01369 {
01370
01371 if ( isProjectile )
01372 {
01373 FirstString= new G4ExcitedString(end,start, +1);
01374 } else {
01375 FirstString= new G4ExcitedString(start,end, -1);
01376 }
01377
01378 SecondString=0;
01379
01380 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01381 FirstString->SetPosition(hadron->GetPosition());
01382
01383
01384
01385 G4double Momentum=hadron->Get4Momentum().vect().mag();
01386 G4double Plus=hadron->Get4Momentum().e() + Momentum;
01387 G4double Minus=hadron->Get4Momentum().e() - Momentum;
01388
01389 G4ThreeVector tmp;
01390 if(Momentum > 0.)
01391 {
01392 tmp.set(hadron->Get4Momentum().px(),
01393 hadron->Get4Momentum().py(),
01394 hadron->Get4Momentum().pz());
01395 tmp/=Momentum;
01396 }
01397 else
01398 {
01399 tmp.set(0.,0.,1.);
01400 }
01401
01402 G4LorentzVector Pstart1(tmp,0.);
01403 G4LorentzVector Pend1(tmp,0.);
01404
01405 if(isProjectile)
01406 {
01407 Pstart1*=(-1.)*Minus/2.;
01408 Pend1 *=(+1.)*Plus /2.;
01409 }
01410 else
01411 {
01412 Pstart1*=(+1.)*Plus/2.;
01413 Pend1 *=(-1.)*Minus/2.;
01414 }
01415
01416 Momentum=-Pstart1.mag();
01417 Pstart1.setT(Momentum);
01418
01419 Momentum=-Pend1.mag();
01420 Pend1.setT(Momentum);
01421
01422 start->Set4Momentum(Pstart1);
01423 end->Set4Momentum(Pend1);
01424 SecondString=0;
01425 }
01426
01427
01428
01429
01430 #ifdef G4_FTFDEBUG
01431 G4cout << " generated string flavors "
01432 << start->GetPDGcode() << " / "
01433 << end->GetPDGcode() << G4endl;
01434 G4cout << " generated string momenta: quark "
01435 << start->Get4Momentum() << "mass : "
01436 <<start->Get4Momentum().mag() << G4endl;
01437 G4cout << " generated string momenta: Diquark "
01438 << end ->Get4Momentum()
01439 << "mass : " <<end->Get4Momentum().mag()<< G4endl;
01440 G4cout << " sum of ends " << Pstart+Pend << G4endl;
01441 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
01442 #endif
01443
01444 return;
01445
01446 }
01447
01448
01449
01450
01451
01452 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
01453 {
01454
01455
01456
01457 G4double range=Pmax-Pmin;
01458
01459 if ( Pmin <= 0. || range <=0. )
01460 {
01461 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
01462 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
01463 }
01464
01465 G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());
01466
01467 return P;
01468 }
01469
01470
01471 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
01472 G4double maxPtSquare) const
01473 {
01474
01475 G4double Pt2(0.);
01476 if(AveragePt2 <= 0.) {Pt2=0.;}
01477 else
01478 {
01479 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
01480 (std::exp(-maxPtSquare/AveragePt2)-1.));
01481 }
01482 G4double Pt=std::sqrt(Pt2);
01483 G4double phi=G4UniformRand() * twopi;
01484 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
01485 }
01486
01487
01488 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
01489 {
01490 G4double z, yf;
01491 do {
01492 z = zmin + G4UniformRand()*(zmax-zmin);
01493 yf = z*z +sqr(1 - z);
01494 }
01495 while (G4UniformRand() > yf);
01496 return z;
01497 }
01498
01499 void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const
01500 {
01501 G4int absIdPDG = std::abs(IdPDG);
01502 Q1= absIdPDG/ 100;
01503 Q2= (absIdPDG %100)/10;
01504
01505 G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
01506
01507 if (IdPDG < 0 ) anti *=-1;
01508 Q1 *= anti;
01509 Q2 *= -1 * anti;
01510 return;
01511 }
01512
01513 void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG,
01514 G4int &Q1, G4int &Q2, G4int &Q3) const
01515 {
01516 Q1 = IdPDG / 1000;
01517 Q2 = (IdPDG % 1000) / 100;
01518 Q3 = (IdPDG % 100) / 10;
01519 return;
01520 }
01521
01522 G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const
01523 {
01524 G4int TmpQ(0);
01525 if( Q3 > Q2 )
01526 {
01527 TmpQ = Q2;
01528 Q2 = Q3;
01529 Q3 = TmpQ;
01530 } else if( Q3 > Q1 )
01531 {
01532 TmpQ = Q1;
01533 Q1 = Q3;
01534 Q3 = TmpQ;
01535 }
01536
01537 if( Q2 > Q1 )
01538 {
01539 TmpQ = Q1;
01540 Q1 = Q2;
01541 Q2 = TmpQ;
01542 }
01543
01544 G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2;
01545 return NewCode;
01546 }
01547
01548 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
01549 {
01550 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
01551 }
01552
01553
01554 G4DiffractiveExcitation::~G4DiffractiveExcitation()
01555 {
01556 }
01557
01558
01559 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
01560 {
01561 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
01562 return *this;
01563 }
01564
01565
01566 int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
01567 {
01568 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
01569 return false;
01570 }
01571
01572 int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
01573 {
01574 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
01575 return true;
01576 }