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 "G4DiffractiveSplitableHadron.hh"
00052 #include "G4DiffractiveExcitation.hh"
00053 #include "G4FTFParameters.hh"
00054 #include "G4ElasticHNScattering.hh"
00055 #include "G4FTFAnnihilation.hh"
00056
00057 #include "G4LorentzRotation.hh"
00058 #include "G4RotationMatrix.hh"
00059 #include "G4ThreeVector.hh"
00060 #include "G4ParticleDefinition.hh"
00061 #include "G4VSplitableHadron.hh"
00062 #include "G4ExcitedString.hh"
00063 #include "G4ParticleTable.hh"
00064 #include "G4Neutron.hh"
00065 #include "G4ParticleDefinition.hh"
00066
00067
00068
00069
00070 G4FTFAnnihilation::G4FTFAnnihilation()
00071 {
00072 }
00073
00074
00075 G4bool G4FTFAnnihilation::
00076 Annihilate(G4VSplitableHadron *projectile,
00077 G4VSplitableHadron *target,
00078 G4VSplitableHadron *&AdditionalString,
00079 G4FTFParameters *theParameters) const
00080 {
00081
00082 G4LorentzVector Pprojectile=projectile->Get4Momentum();
00083
00084
00085
00086 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
00087 if(ProjectilePDGcode > 0)
00088 {
00089 target->SetStatus(2);
00090 return false;
00091 }
00092
00093
00094
00095
00096 G4double M0projectile2=Pprojectile.mag2();
00097
00098 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
00099
00100 G4LorentzVector Ptarget=target->Get4Momentum();
00101
00102
00103
00104
00105
00106
00107 G4double M0target2=Ptarget.mag2();
00108
00109
00110
00111
00112
00113
00114 G4double AveragePt2=theParameters->GetAveragePt2();
00115
00116
00117 G4LorentzVector Psum;
00118 Psum=Pprojectile+Ptarget;
00119 G4double S=Psum.mag2();
00120
00121
00122 G4LorentzRotation toCms(-1*Psum.boostVector());
00123
00124 G4LorentzVector Ptmp=toCms*Pprojectile;
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135 toCms.rotateZ(-1*Ptmp.phi());
00136 toCms.rotateY(-1*Ptmp.theta());
00137
00138 G4LorentzRotation toLab(toCms.inverse());
00139
00140 G4double SqrtS=std::sqrt(S);
00141
00142 G4double maxPtSquare;
00143
00144
00145
00146 G4double X_a(0.), X_b(0.), X_c(0.), X_d(0.);
00147 G4double MesonProdThreshold=projectile->GetDefinition()->GetPDGMass()+
00148 target->GetDefinition()->GetPDGMass()+
00149 (2.*140.+16.)*MeV;
00150
00151 G4double Prel2= S*S + M0projectile2*M0projectile2 + M0target2*M0target2 -
00152 2.*S*M0projectile2 - 2.*S*M0target2 - 2.*M0projectile2*M0target2;
00153 Prel2/=S;
00154
00155 if(Prel2 < 0. )
00156 {
00157 X_a= 625.1;
00158 X_b= 9.780;
00159 X_c= 49.989;
00160 X_d= 6.614;
00161 }
00162 else
00163 {
00164 G4double FlowF=1./std::sqrt(Prel2)*GeV;
00165
00166
00167
00168
00169 X_a=25.*FlowF;
00170
00171
00172 if(SqrtS < MesonProdThreshold)
00173 {
00174 X_b=3.13+140.*std::pow((MesonProdThreshold - SqrtS)/GeV,2.5);
00175 }
00176 else
00177 {
00178 X_b=6.8*GeV/SqrtS;
00179 }
00180 if(projectile->GetDefinition()->GetPDGMass()+
00181 target->GetDefinition()->GetPDGMass() > SqrtS) {X_b=0.;}
00182
00183
00184
00185 X_c=2.*FlowF*sqr(projectile->GetDefinition()->GetPDGMass()+
00186 target->GetDefinition()->GetPDGMass())/S;
00187
00188
00189 X_d=23.3*GeV*GeV/S;
00190 }
00191
00192
00193
00194 if((ProjectilePDGcode == -2212)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00195 {X_b*=5.; X_c*=5.; X_d*=6.;}
00196 else if((ProjectilePDGcode == -2212)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00197 {X_b*=4.; X_c*=4.; X_d*=4.;}
00198 else if((ProjectilePDGcode == -2112)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00199 {X_b*=4.; X_c*=4.; X_d*=4.;}
00200 else if((ProjectilePDGcode == -2112)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00201 {X_b*=5.; X_c*=5.; X_d*=6.;}
00202 else if((ProjectilePDGcode == -3122)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00203 {X_b*=3.; X_c*=3.; X_d*=2.;}
00204 else if((ProjectilePDGcode == -3122)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00205 {X_b*=3.; X_c*=3.; X_d*=2.;}
00206 else if((ProjectilePDGcode == -3112)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00207 {X_b*=2.; X_c*=2.; X_d*=0.;}
00208 else if((ProjectilePDGcode == -3112)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00209 {X_b*=4.; X_c*=4.; X_d*=2.;}
00210 else if((ProjectilePDGcode == -3212)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00211 {X_b*=3.; X_c*=3.; X_d*=2.;}
00212 else if((ProjectilePDGcode == -3212)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00213 {X_b*=3.; X_c*=3.; X_d*=2.;}
00214 else if((ProjectilePDGcode == -3222)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00215 {X_b*=4.; X_c*=4.; X_d*=2.;}
00216 else if((ProjectilePDGcode == -3222)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00217 {X_b*=2.; X_c*=2.; X_d*=0.;}
00218 else if((ProjectilePDGcode == -3312)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00219 {X_b*=1.; X_c*=1.; X_d*=0.;}
00220 else if((ProjectilePDGcode == -3312)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00221 {X_b*=2.; X_c*=2.; X_d*=0.;}
00222 else if((ProjectilePDGcode == -3322)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00223 {X_b*=2.; X_c*=2.; X_d*=0.;}
00224 else if((ProjectilePDGcode == -3322)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00225 {X_b*=1.; X_c*=1.; X_d*=0.;}
00226 else if((ProjectilePDGcode == -3334)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
00227 {X_b*=0.; X_c*=0.; X_d*=0.;}
00228 else if((ProjectilePDGcode == -3334)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
00229 {X_b*=0.; X_c*=0.; X_d*=0.;}
00230 else {G4cout<<"Unknown anti-baryon for FTF annihilation: PDGcodes - "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;}
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 G4double Xannihilation=X_a+X_b+X_c+X_d;
00241
00242
00243 G4int AQ[3];
00244 UnpackBaryon(ProjectilePDGcode, AQ[0], AQ[1], AQ[2]);
00245
00246
00247 G4int Q[3];
00248 UnpackBaryon(TargetPDGcode, Q[0], Q[1], Q[2]);
00249
00250
00251 G4double Ksi=G4UniformRand();
00252
00253 if(Ksi < X_a/Xannihilation)
00254 {
00255
00256
00257
00258
00259 G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(6));
00260
00261 G4int Tmp1(0), Tmp2(0);
00262 if(SampledCase == 0) { }
00263 if(SampledCase == 1) {Tmp1=AQ[1]; AQ[1]=AQ[2]; AQ[2]=Tmp1;}
00264 if(SampledCase == 2) {Tmp1=AQ[0]; AQ[0]=AQ[1]; AQ[1]=Tmp1;}
00265 if(SampledCase == 3) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=AQ[2]; AQ[1]=Tmp1; AQ[2]=Tmp2;}
00266 if(SampledCase == 4) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=Tmp2; AQ[1]=AQ[2]; AQ[2]=Tmp1;}
00267 if(SampledCase == 5) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=AQ[2]; AQ[1]=Tmp2; AQ[2]=Tmp1;}
00268
00269
00270
00271 projectile->SplitUp();
00272
00273 projectile->SetFirstParton(AQ[0]);
00274 projectile->SetSecondParton(Q[0]);
00275 projectile->SetStatus(1);
00276
00277
00278 target->SplitUp();
00279
00280 target->SetFirstParton(Q[1]);
00281 target->SetSecondParton(AQ[1]);
00282 target->SetStatus(1);
00283
00284
00285 AdditionalString=new G4DiffractiveSplitableHadron();
00286 AdditionalString->SplitUp();
00287 AdditionalString->SetFirstParton(AQ[2]);
00288 AdditionalString->SetSecondParton(Q[2]);
00289 AdditionalString->SetStatus(1);
00290
00291
00292
00293
00294
00295 G4ThreeVector Quark_Mom[6];
00296 G4double ModMom2[6];
00297
00298
00299 AveragePt2=200.*200.; maxPtSquare=S;
00300
00301
00302 G4double SumMt(0.);
00303
00304 G4double MassQ2=0.;
00305
00306 G4int NumberOfTries(0);
00307 G4double ScaleFactor(1.);
00308 do
00309 {
00310 NumberOfTries++;
00311
00312 if(NumberOfTries == 100*(NumberOfTries/100))
00313 {
00314 ScaleFactor/=2.;
00315 AveragePt2 *=ScaleFactor;
00316 }
00317
00318 G4ThreeVector PtSum(0.,0.,0.);
00319 for(G4int i=0; i<6; i++)
00320 {
00321 Quark_Mom[i]=GaussianPt(AveragePt2, maxPtSquare);
00322 PtSum+=Quark_Mom[i];
00323 }
00324
00325 PtSum/=6.;
00326
00327 SumMt=0.;
00328 for(G4int i=0; i<6; i++)
00329 {
00330 Quark_Mom[i]-=PtSum;
00331
00332 ModMom2[i]=Quark_Mom[i].mag2();
00333 SumMt+=std::sqrt(ModMom2[i]+MassQ2);
00334 }
00335 } while(SumMt > SqrtS);
00336
00337 G4double WminusTarget(0.), WplusProjectile(0.);
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 G4double Alfa_R=0.5;
00362
00363 NumberOfTries=0;
00364 ScaleFactor=1.;
00365
00366 G4bool Succes(true);
00367 do
00368 {
00369 Succes=true;
00370 NumberOfTries++;
00371
00372 if(NumberOfTries == 100*(NumberOfTries/100))
00373 {
00374 ScaleFactor/=2.;
00375 }
00376
00377 if(Alfa_R == 1.)
00378 {
00379 G4double Xaq1=1.-std::sqrt(G4UniformRand());
00380 G4double Xaq2=(1.-Xaq1)*G4UniformRand();
00381 G4double Xaq3=1.-Xaq1-Xaq2;
00382
00383 Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2); Quark_Mom[2].setZ(Xaq3);
00384 }
00385 else
00386 {
00387 G4double Xaq1=sqr(G4UniformRand());
00388 G4double Xaq2=(1.-Xaq1)*sqr(std::sin(pi/2.*G4UniformRand()));
00389 G4double Xaq3=1.-Xaq1-Xaq2;
00390
00391 Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2); Quark_Mom[2].setZ(Xaq3);
00392 }
00393
00394
00395 if(Alfa_R == 1.)
00396 {
00397 G4double Xq1=1.-std::sqrt(G4UniformRand());
00398 G4double Xq2=(1.-Xq1)*G4UniformRand();
00399 G4double Xq3=1.-Xq1-Xq2;
00400
00401 Quark_Mom[3].setZ(Xq1); Quark_Mom[4].setZ(Xq2); Quark_Mom[5].setZ(Xq3);
00402 }
00403 else
00404 {
00405 G4double Xq1=sqr(G4UniformRand());
00406 G4double Xq2=(1.-Xq1)*sqr(std::sin(pi/2.*G4UniformRand()));
00407 G4double Xq3=1.-Xq1-Xq2;
00408
00409 Quark_Mom[3].setZ(Xq1); Quark_Mom[4].setZ(Xq2); Quark_Mom[5].setZ(Xq3);
00410 }
00411
00412 G4double Alfa(0.), Beta(0.);
00413
00414 for(G4int i=0; i<3; i++)
00415 {
00416 if(Quark_Mom[i].getZ() != 0.)
00417 {Alfa+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
00418 else {Succes=false;}
00419 }
00420
00421 for(G4int i=3; i<6; i++)
00422 {
00423 if(Quark_Mom[i].getZ() != 0.)
00424 {Beta+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
00425 else {Succes=false;}
00426 }
00427
00428 if(!Succes) continue;
00429
00430 if(std::sqrt(Alfa)+std::sqrt(Beta) > SqrtS) {Succes=false; continue;}
00431
00432 G4double DecayMomentum2=S*S + Alfa*Alfa + Beta*Beta
00433 - 2.*S*Alfa - 2.*S*Beta - 2.*Alfa*Beta;
00434
00435 WminusTarget=(S-Alfa+Beta+std::sqrt(DecayMomentum2))/2./SqrtS;
00436 WplusProjectile=SqrtS-Beta/WminusTarget;
00437
00438 } while(!Succes);
00439
00440
00441 G4double SqrtScaleF=std::sqrt(ScaleFactor);
00442
00443 for(G4int i=0; i<3; i++)
00444 {
00445 G4double Pz=WplusProjectile*Quark_Mom[i].getZ()/2.-
00446 (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WplusProjectile*Quark_Mom[i].getZ());
00447 Quark_Mom[i].setZ(Pz);
00448
00449 if(ScaleFactor != 1.)
00450 {
00451 Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
00452 Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
00453 }
00454 }
00455
00456 for(G4int i=3; i<6; i++)
00457 {
00458 G4double Pz=-WminusTarget*Quark_Mom[i].getZ()/2.+
00459 (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WminusTarget*Quark_Mom[i].getZ());
00460 Quark_Mom[i].setZ(Pz);
00461
00462 if(ScaleFactor != 1.)
00463 {
00464 Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
00465 Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
00466 }
00467 }
00468
00469
00470
00471
00472 G4ThreeVector tmp=Quark_Mom[0]+Quark_Mom[3];
00473 G4LorentzVector Pstring1(tmp,std::sqrt(Quark_Mom[0].mag2()+MassQ2)+
00474 std::sqrt(Quark_Mom[3].mag2()+MassQ2));
00475 G4double Ystring1=Pstring1.rapidity();
00476
00477
00478
00479
00480
00481
00482
00483
00484 tmp=Quark_Mom[1]+Quark_Mom[4];
00485 G4LorentzVector Pstring2(tmp,std::sqrt(Quark_Mom[1].mag2()+MassQ2)+
00486 std::sqrt(Quark_Mom[4].mag2()+MassQ2));
00487 G4double Ystring2=Pstring2.rapidity();
00488
00489
00490
00491
00492
00493
00494
00495
00496 tmp=Quark_Mom[2]+Quark_Mom[5];
00497 G4LorentzVector Pstring3(tmp,std::sqrt(Quark_Mom[2].mag2()+MassQ2)+
00498 std::sqrt(Quark_Mom[5].mag2()+MassQ2));
00499 G4double Ystring3=Pstring3.rapidity();
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 G4LorentzVector LeftString(0.,0.,0.,0.);
00512
00513 if((Ystring1 > Ystring2)&&(Ystring2 > Ystring3))
00514 {
00515 Pprojectile=Pstring1;
00516 LeftString =Pstring2;
00517 Ptarget =Pstring3;
00518 }
00519
00520 if((Ystring1 > Ystring3)&&(Ystring3 > Ystring2))
00521 {
00522 Pprojectile=Pstring1;
00523 LeftString =Pstring3;
00524 Ptarget =Pstring2;
00525 }
00526
00527 if((Ystring2 > Ystring1)&&(Ystring1 > Ystring3))
00528 {
00529 Pprojectile=Pstring2;
00530 LeftString =Pstring1;
00531 Ptarget =Pstring3;
00532 }
00533
00534 if((Ystring2 > Ystring3)&&(Ystring3 > Ystring1))
00535 {
00536 Pprojectile=Pstring2;
00537 LeftString =Pstring3;
00538 Ptarget =Pstring1;
00539 }
00540
00541 if((Ystring3 > Ystring1)&&(Ystring1 > Ystring2))
00542 {
00543 Pprojectile=Pstring3;
00544 LeftString =Pstring1;
00545 Ptarget =Pstring2;
00546 }
00547
00548 if((Ystring3 > Ystring2)&&(Ystring2 > Ystring1))
00549 {
00550 Pprojectile=Pstring3;
00551 LeftString =Pstring2;
00552 Ptarget =Pstring1;
00553 }
00554
00555
00556
00557
00558 Pprojectile.transform(toLab);
00559 LeftString.transform(toLab);
00560 Ptarget.transform(toLab);
00561
00562
00563
00564 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00565 projectile->SetPosition(target->GetPosition());
00566
00567 AdditionalString->SetTimeOfCreation(target->GetTimeOfCreation());
00568 AdditionalString->SetPosition(target->GetPosition());
00569
00570
00571
00572
00573
00574
00575 projectile->Set4Momentum(Pprojectile);
00576 AdditionalString->Set4Momentum(LeftString);
00577 target->Set4Momentum(Ptarget);
00578
00579 projectile->IncrementCollisionCount(1);
00580 AdditionalString->IncrementCollisionCount(1);
00581 target->IncrementCollisionCount(1);
00582
00583 return true;
00584 }
00585
00586
00587
00588
00589 if(Ksi < (X_a+X_b)/Xannihilation)
00590 {
00591
00592 G4int CandidatsN(0), CandAQ[9][2], CandQ[9][2];
00593 G4int LeftAQ1(0), LeftAQ2(0), LeftQ1(0), LeftQ2(0);
00594
00595 for(G4int iAQ=0; iAQ<3; iAQ++)
00596 {
00597 for(G4int iQ=0; iQ<3; iQ++)
00598 {
00599 if(-AQ[iAQ] == Q[iQ])
00600 {
00601 if(iAQ == 0) {CandAQ[CandidatsN][0]=1; CandAQ[CandidatsN][1]=2;}
00602 if(iAQ == 1) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=2;}
00603 if(iAQ == 2) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=1;}
00604 if(iQ == 0) {CandQ[CandidatsN][0] =1; CandQ[CandidatsN][1]=2;}
00605 if(iQ == 1) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=2;}
00606 if(iQ == 2) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=1;}
00607 CandidatsN++;
00608 }
00609 }
00610 }
00611
00612
00613
00614 if(CandidatsN != 0)
00615 {
00616 G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
00617
00618 LeftAQ1=AQ[CandAQ[SampledCase][0]];
00619 LeftAQ2=AQ[CandAQ[SampledCase][1]];
00620
00621 LeftQ1=Q[CandQ[SampledCase][0]];
00622 LeftQ2=Q[CandQ[SampledCase][1]];
00623
00624
00625 G4int Anti_DQ(0), DQ(0);
00626
00627 if(std::abs(LeftAQ1) > std::abs(LeftAQ2))
00628 {
00629 Anti_DQ=1000*LeftAQ1+100*LeftAQ2-3;
00630 } else
00631 {
00632 Anti_DQ=1000*LeftAQ2+100*LeftAQ1-3;
00633 }
00634
00635
00636 if(std::abs(LeftQ1) > std::abs(LeftQ2))
00637 {
00638 DQ=1000*LeftQ1+100*LeftQ2+3;
00639 } else
00640 {
00641 DQ=1000*LeftQ2+100*LeftQ1+3;
00642 }
00643
00644
00645
00646
00647
00648 projectile->SplitUp();
00649
00650
00651
00652 projectile->SetFirstParton(DQ);
00653 projectile->SetSecondParton(Anti_DQ);
00654
00655 projectile->SetStatus(1);
00656 target->SetStatus(3);
00657
00658 Pprojectile.setPx(0.);
00659 Pprojectile.setPy(0.);
00660 Pprojectile.setPz(0.);
00661 Pprojectile.setE(SqrtS);
00662 Pprojectile.transform(toLab);
00663
00664
00665 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00666 projectile->SetPosition(target->GetPosition());
00667
00668
00669
00670
00671
00672
00673 projectile->Set4Momentum(Pprojectile);
00674
00675 projectile->IncrementCollisionCount(1);
00676
00677 return true;
00678 }
00679 }
00680
00681
00682
00683 if(Ksi < (X_a+X_b+X_c)/Xannihilation)
00684 {
00685
00686
00687
00688 G4int CandidatsN(0), CandAQ[9][2], CandQ[9][2];
00689 G4int LeftAQ1(0), LeftAQ2(0), LeftQ1(0), LeftQ2(0);
00690
00691 for(G4int iAQ=0; iAQ<3; iAQ++)
00692 {
00693 for(G4int iQ=0; iQ<3; iQ++)
00694 {
00695 if(-AQ[iAQ] == Q[iQ])
00696 {
00697 if(iAQ == 0) {CandAQ[CandidatsN][0]=1; CandAQ[CandidatsN][1]=2;}
00698 if(iAQ == 1) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=2;}
00699 if(iAQ == 2) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=1;}
00700 if(iQ == 0) {CandQ[CandidatsN][0] =1; CandQ[CandidatsN][1]=2;}
00701 if(iQ == 1) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=2;}
00702 if(iQ == 2) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=1;}
00703 CandidatsN++;
00704 }
00705 }
00706 }
00707
00708
00709
00710 if(CandidatsN != 0)
00711 {
00712 G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
00713
00714 LeftAQ1=AQ[CandAQ[SampledCase][0]];
00715 LeftAQ2=AQ[CandAQ[SampledCase][1]];
00716
00717 if(G4UniformRand() < 0.5)
00718 {
00719 LeftQ1=Q[CandQ[SampledCase][0]];
00720 LeftQ2=Q[CandQ[SampledCase][1]];
00721 } else
00722 {
00723 LeftQ2=Q[CandQ[SampledCase][0]];
00724 LeftQ1=Q[CandQ[SampledCase][1]];
00725 }
00726
00727
00728
00729 projectile->SplitUp();
00730
00731 projectile->SetFirstParton(LeftAQ1);
00732 projectile->SetSecondParton(LeftQ1);
00733 projectile->SetStatus(1);
00734
00735
00736 target->SplitUp();
00737
00738 target->SetFirstParton(LeftQ2);
00739 target->SetSecondParton(LeftAQ2);
00740 target->SetStatus(1);
00741
00742
00743
00744
00745 G4ThreeVector Quark_Mom[4];
00746 G4double ModMom2[4];
00747
00748
00749 AveragePt2=200.*200.; maxPtSquare=S;
00750
00751
00752 G4double SumMt(0.);
00753
00754 G4double MassQ2=0.;
00755
00756 G4int NumberOfTries(0);
00757 G4double ScaleFactor(1.);
00758 do
00759 {
00760 NumberOfTries++;
00761
00762 if(NumberOfTries == 100*(NumberOfTries/100))
00763 {
00764 ScaleFactor/=2.;
00765 AveragePt2 *=ScaleFactor;
00766 }
00767
00768 G4ThreeVector PtSum(0.,0.,0.);
00769 for(G4int i=0; i<4; i++)
00770 {
00771 Quark_Mom[i]=GaussianPt(AveragePt2, maxPtSquare);
00772 PtSum+=Quark_Mom[i];
00773 }
00774
00775 PtSum/=4.;
00776
00777 SumMt=0.;
00778 for(G4int i=0; i<4; i++)
00779 {
00780 Quark_Mom[i]-=PtSum;
00781
00782 ModMom2[i]=Quark_Mom[i].mag2();
00783 SumMt+=std::sqrt(ModMom2[i]+MassQ2);
00784 }
00785 } while(SumMt > SqrtS);
00786
00787 G4double WminusTarget(0.), WplusProjectile(0.);
00788
00789
00790 G4double Alfa_R=0.5;
00791
00792 NumberOfTries=0;
00793 ScaleFactor=1.;
00794
00795 G4bool Succes(true);
00796 do
00797 {
00798 Succes=true;
00799 NumberOfTries++;
00800
00801 if(NumberOfTries == 100*(NumberOfTries/100))
00802 {
00803 ScaleFactor/=2.;
00804 }
00805
00806 if(Alfa_R == 1.)
00807 {
00808 G4double Xaq1=std::sqrt(G4UniformRand());
00809 G4double Xaq2=1.-Xaq1;
00810
00811 Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2);
00812 }
00813 else
00814 {
00815 G4double Xaq1=sqr(std::sin(pi/2.*G4UniformRand()));
00816 G4double Xaq2=1.-Xaq1;
00817
00818 Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2);
00819 }
00820
00821
00822 if(Alfa_R == 1.)
00823 {
00824 G4double Xq1=1.-std::sqrt(G4UniformRand());
00825 G4double Xq2=1.-Xq1;
00826
00827 Quark_Mom[2].setZ(Xq1); Quark_Mom[3].setZ(Xq2);
00828 }
00829 else
00830 {
00831 G4double Xq1=sqr(std::sin(pi/2.*G4UniformRand()));
00832 G4double Xq2=1.-Xq1;
00833
00834 Quark_Mom[2].setZ(Xq1); Quark_Mom[3].setZ(Xq2);
00835 }
00836
00837 G4double Alfa(0.), Beta(0.);
00838
00839 for(G4int i=0; i<2; i++)
00840 {
00841 if(Quark_Mom[i].getZ() != 0.)
00842 {Alfa+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
00843 else {Succes=false;}
00844 }
00845
00846 for(G4int i=2; i<4; i++)
00847 {
00848 if(Quark_Mom[i].getZ() != 0.)
00849 {Beta+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
00850 else {Succes=false;}
00851 }
00852
00853 if(!Succes) continue;
00854
00855 if(std::sqrt(Alfa)+std::sqrt(Beta) > SqrtS) {Succes=false; continue;}
00856
00857 G4double DecayMomentum2=S*S + Alfa*Alfa + Beta*Beta
00858 - 2.*S*Alfa - 2.*S*Beta - 2.*Alfa*Beta;
00859
00860 WminusTarget=(S-Alfa+Beta+std::sqrt(DecayMomentum2))/2./SqrtS;
00861 WplusProjectile=SqrtS-Beta/WminusTarget;
00862
00863 } while(!Succes);
00864
00865
00866 G4double SqrtScaleF=std::sqrt(ScaleFactor);
00867
00868 for(G4int i=0; i<2; i++)
00869 {
00870 G4double Pz=WplusProjectile*Quark_Mom[i].getZ()/2.-
00871 (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WplusProjectile*Quark_Mom[i].getZ());
00872 Quark_Mom[i].setZ(Pz);
00873
00874 if(ScaleFactor != 1.)
00875 {
00876 Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
00877 Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
00878 }
00879
00880 }
00881
00882 for(G4int i=2; i<4; i++)
00883 {
00884 G4double Pz=-WminusTarget*Quark_Mom[i].getZ()/2.+
00885 (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WminusTarget*Quark_Mom[i].getZ());
00886 Quark_Mom[i].setZ(Pz);
00887
00888 if(ScaleFactor != 1.)
00889 {
00890 Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
00891 Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
00892 }
00893
00894 }
00895
00896
00897
00898
00899 G4ThreeVector tmp=Quark_Mom[0]+Quark_Mom[2];
00900 G4LorentzVector Pstring1(tmp,std::sqrt(Quark_Mom[0].mag2()+MassQ2)+
00901 std::sqrt(Quark_Mom[2].mag2()+MassQ2));
00902 G4double Ystring1=Pstring1.rapidity();
00903
00904
00905
00906
00907
00908
00909
00910
00911 tmp=Quark_Mom[1]+Quark_Mom[3];
00912 G4LorentzVector Pstring2(tmp,std::sqrt(Quark_Mom[1].mag2()+MassQ2)+
00913 std::sqrt(Quark_Mom[3].mag2()+MassQ2));
00914 G4double Ystring2=Pstring2.rapidity();
00915
00916
00917
00918
00919
00920
00921
00922
00923 if(Ystring1 > Ystring2)
00924 {
00925 Pprojectile=Pstring1;
00926 Ptarget =Pstring2;
00927 } else
00928 {
00929 Pprojectile=Pstring2;
00930 Ptarget =Pstring1;
00931 }
00932
00933
00934
00935
00936 Pprojectile.transform(toLab);
00937 Ptarget.transform(toLab);
00938
00939
00940
00941 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00942 projectile->SetPosition(target->GetPosition());
00943
00944
00945
00946
00947
00948
00949
00950 projectile->Set4Momentum(Pprojectile);
00951
00952 target->Set4Momentum(Ptarget);
00953
00954 projectile->IncrementCollisionCount(1);
00955 target->IncrementCollisionCount(1);
00956
00957 return true;
00958 }
00959 }
00960
00961
00962
00963
00964 if(Ksi < (X_a+X_b+X_c+X_d)/Xannihilation)
00965 {
00966
00967 G4int CandidatsN(0), CandAQ[9], CandQ[9];
00968 G4int LeftAQ(0), LeftQ(0);
00969
00970 for(G4int iAQ1=0; iAQ1<3; iAQ1++)
00971 {
00972 for(G4int iAQ2=0; iAQ2<3; iAQ2++)
00973 {
00974 if(iAQ1 != iAQ2)
00975 {
00976 for(G4int iQ1=0; iQ1<3; iQ1++)
00977 {
00978 for(G4int iQ2=0; iQ2<3; iQ2++)
00979 {
00980 if(iQ1 != iQ2)
00981 {
00982 if((-AQ[iAQ1] == Q[iQ1]) && (-AQ[iAQ2] == Q[iQ2]))
00983 {
00984 if((iAQ1 == 0) && (iAQ2 == 1)){CandAQ[CandidatsN]=2;}
00985 if((iAQ1 == 1) && (iAQ2 == 0)){CandAQ[CandidatsN]=2;}
00986
00987 if((iAQ1 == 0) && (iAQ2 == 2)){CandAQ[CandidatsN]=1;}
00988 if((iAQ1 == 2) && (iAQ2 == 0)){CandAQ[CandidatsN]=1;}
00989
00990 if((iAQ1 == 1) && (iAQ2 == 2)){CandAQ[CandidatsN]=0;}
00991 if((iAQ1 == 2) && (iAQ2 == 1)){CandAQ[CandidatsN]=0;}
00992
00993 if((iQ1 == 0) && (iQ2 == 1)){CandQ[CandidatsN]=2;}
00994 if((iQ1 == 1) && (iQ2 == 0)){CandQ[CandidatsN]=2;}
00995
00996 if((iQ1 == 0) && (iQ2 == 2)){CandQ[CandidatsN]=1;}
00997 if((iQ1 == 2) && (iQ2 == 0)){CandQ[CandidatsN]=1;}
00998
00999 if((iQ1 == 1) && (iQ2 == 2)){CandQ[CandidatsN]=0;}
01000 if((iQ1 == 2) && (iQ2 == 1)){CandQ[CandidatsN]=0;}
01001 CandidatsN++;
01002 }
01003 }
01004 }
01005 }
01006 }
01007 }
01008 }
01009
01010
01011 if(CandidatsN != 0)
01012 {
01013 G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
01014
01015 LeftAQ=AQ[CandAQ[SampledCase]];
01016
01017 LeftQ =Q[CandQ[SampledCase]];
01018
01019
01020
01021
01022 projectile->SplitUp();
01023
01024
01025
01026 projectile->SetFirstParton(LeftQ);
01027 projectile->SetSecondParton(LeftAQ);
01028
01029 projectile->SetStatus(1);
01030 target->SetStatus(3);
01031
01032 Pprojectile.setPx(0.);
01033 Pprojectile.setPy(0.);
01034 Pprojectile.setPz(0.);
01035 Pprojectile.setE(SqrtS);
01036 Pprojectile.transform(toLab);
01037
01038
01039 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
01040 projectile->SetPosition(target->GetPosition());
01041
01042
01043
01044
01045
01046
01047 projectile->Set4Momentum(Pprojectile);
01048
01049 projectile->IncrementCollisionCount(1);
01050 return true;
01051 }
01052 }
01053
01054
01055
01056 return true;
01057 }
01058
01059
01060 G4double G4FTFAnnihilation::ChooseX(G4double Alpha, G4double Beta) const
01061 {
01062
01063
01064 G4double tmp=Alpha*Beta;
01065 tmp*=1.;
01066 return 0.5;
01067 }
01068
01069
01070 G4ThreeVector G4FTFAnnihilation::GaussianPt(G4double AveragePt2,
01071 G4double maxPtSquare) const
01072 {
01073
01074 G4double Pt2(0.);
01075 if(AveragePt2 <= 0.) {Pt2=0.;}
01076 else
01077 {
01078 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
01079 (std::exp(-maxPtSquare/AveragePt2)-1.));
01080 }
01081 G4double Pt=std::sqrt(Pt2);
01082 G4double phi=G4UniformRand() * twopi;
01083 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
01084 }
01085
01086
01087
01088 void G4FTFAnnihilation::UnpackBaryon(G4int IdPDG,
01089 G4int &Q1, G4int &Q2, G4int &Q3) const
01090 {
01091 G4int AbsId=std::abs(IdPDG);
01092
01093 Q1 = AbsId / 1000;
01094 Q2 = (AbsId % 1000) / 100;
01095 Q3 = (AbsId % 100) / 10;
01096
01097 if(IdPDG < 0 ) {Q1=-Q1; Q2=-Q2; Q3=-Q3;}
01098
01099 return;
01100 }
01101
01102
01103 G4FTFAnnihilation::G4FTFAnnihilation(const G4FTFAnnihilation &)
01104 {
01105 throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation copy contructor not meant to be called");
01106 }
01107
01108
01109 G4FTFAnnihilation::~G4FTFAnnihilation()
01110 {
01111 }
01112
01113
01114 const G4FTFAnnihilation & G4FTFAnnihilation::operator=(const G4FTFAnnihilation &)
01115 {
01116 throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation = operator not meant to be called");
01117
01118 }
01119
01120
01121 int G4FTFAnnihilation::operator==(const G4FTFAnnihilation &) const
01122 {
01123 throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation == operator not meant to be called");
01124
01125 }
01126
01127 int G4FTFAnnihilation::operator!=(const G4FTFAnnihilation &) const
01128 {
01129 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator not meant to be called");
01130
01131 }