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 #include "G4HEInelastic.hh"
00043 #include "globals.hh"
00044 #include "G4ios.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4HEVector.hh"
00048 #include "G4ParticleDefinition.hh"
00049 #include "G4DynamicParticle.hh"
00050 #include "G4ParticleTable.hh"
00051 #include "G4KaonZero.hh"
00052 #include "G4AntiKaonZero.hh"
00053 #include "G4Deuteron.hh"
00054 #include "G4Triton.hh"
00055 #include "G4Alpha.hh"
00056
00057 void G4HEInelastic::FillParticleChange(G4HEVector pv[], G4int aVecLength)
00058 {
00059 theParticleChange.Clear();
00060 for (G4int i=0; i<aVecLength; i++)
00061 {
00062 G4int pdgCode = pv[i].getCode();
00063 G4ParticleDefinition * aDefinition=NULL;
00064 if(pdgCode == 0)
00065 {
00066 G4int bNumber = pv[i].getBaryonNumber();
00067 if(bNumber==2) aDefinition = G4Deuteron::Deuteron();
00068 if(bNumber==3) aDefinition = G4Triton::Triton();
00069 if(bNumber==4) aDefinition = G4Alpha::Alpha();
00070 }
00071 else
00072 {
00073 aDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
00074 }
00075 G4DynamicParticle * aParticle = new G4DynamicParticle();
00076 aParticle->SetDefinition(aDefinition);
00077 aParticle->SetMomentum(pv[i].getMomentum()*GeV);
00078 theParticleChange.AddSecondary(aParticle);
00079 }
00080 }
00081
00082 void G4HEInelastic::SetParticles()
00083 {
00084 PionPlus.setDefinition("PionPlus");
00085 PionZero.setDefinition("PionZero");
00086 PionMinus.setDefinition("PionMinus");
00087 KaonPlus.setDefinition("KaonPlus");
00088 KaonZero.setDefinition("KaonZero");
00089 AntiKaonZero.setDefinition("AntiKaonZero");
00090 KaonMinus.setDefinition("KaonMinus");
00091 KaonZeroShort.setDefinition("KaonZeroShort");
00092 KaonZeroLong.setDefinition("KaonZeroLong");
00093 Proton.setDefinition("Proton");
00094 AntiProton.setDefinition("AntiProton");
00095 Neutron.setDefinition("Neutron");
00096 AntiNeutron.setDefinition("AntiNeutron");
00097 Lambda.setDefinition("Lambda");
00098 AntiLambda.setDefinition("AntiLambda");
00099 SigmaPlus.setDefinition("SigmaPlus");
00100 SigmaZero.setDefinition("SigmaZero");
00101 SigmaMinus.setDefinition("SigmaMinus");
00102 AntiSigmaPlus.setDefinition("AntiSigmaPlus");
00103 AntiSigmaZero.setDefinition("AntiSigmaZero");
00104 AntiSigmaMinus.setDefinition("AntiSigmaMinus");
00105 XiZero.setDefinition("XiZero");
00106 XiMinus.setDefinition("XiMinus");
00107 AntiXiZero.setDefinition("AntiXiZero");
00108 AntiXiMinus.setDefinition("AntiXiMinus");
00109 OmegaMinus.setDefinition("OmegaMinus");
00110 AntiOmegaMinus.setDefinition("AntiOmegaMinus");
00111 Deuteron.setDefinition("Deuteron");
00112 Triton.setDefinition("Triton");
00113 Alpha.setDefinition("Alpha");
00114 Gamma.setDefinition("Gamma");
00115 return;
00116 }
00117
00118 G4double G4HEInelastic::Amin(G4double a, G4double b)
00119 {
00120 G4double c = a;
00121 if(b < a) c = b;
00122 return c;
00123 }
00124
00125 G4double G4HEInelastic::Amax(G4double a, G4double b)
00126 {
00127 G4double c = a;
00128 if(b > a) c = b;
00129 return c;
00130 }
00131
00132 G4int
00133 G4HEInelastic::Imin(G4int a, G4int b)
00134 {
00135 G4int c = a;
00136 if(b < a) c = b;
00137 return c;
00138 }
00139 G4int
00140 G4HEInelastic::Imax(G4int a, G4int b)
00141 {
00142 G4int c = a;
00143 if(b > a) c = b;
00144 return c;
00145 }
00146
00147
00148 G4double
00149 G4HEInelastic::NuclearInelasticity(G4double incidentKineticEnergy,
00150 G4double atomicWeight,
00151 G4double )
00152 {
00153 G4double expu = 82.;
00154 G4double expl = -expu;
00155 G4double ala = std::log(atomicWeight);
00156 G4double ale = std::log(incidentKineticEnergy);
00157 G4double sig1 = 0.5;
00158 G4double sig2 = 0.5;
00159 G4double em = Amin(0.239 + 0.0408*ala*ala, 1.);
00160 G4double cinem = Amin(0.0019*std::pow(ala,3.), 0.15);
00161 G4double sig = (ale > em) ? sig2 : sig1;
00162 G4double corr = Amin(Amax(-std::pow(ale-em,2.)/(2.*sig*sig),expl), expu);
00163 G4double dum1 = -(incidentKineticEnergy)*cinem;
00164 G4double dum2 = std::abs(dum1);
00165 G4double dum3 = std::exp(corr);
00166 G4double cinema = 0.;
00167 if (dum2 >= 1.) cinema = dum1*dum3;
00168 else if (dum3 > 1.e-10) cinema = dum1*dum3;
00169 cinema = - Amax(-incidentKineticEnergy, cinema);
00170 if(verboseLevel > 1) {
00171 G4cout << " NuclearInelasticity: " << ala << " " << ale << " "
00172 << em << " " << corr << " " << dum1 << " " << dum2 << " "
00173 << dum3 << " " << cinema << G4endl;
00174 }
00175 return cinema;
00176 }
00177
00178 G4double
00179 G4HEInelastic::NuclearExcitation(G4double incidentKineticEnergy,
00180 G4double atomicWeight,
00181 G4double atomicNumber,
00182 G4double& excitationEnergyGPN,
00183 G4double& excitationEnergyDTA)
00184 {
00185 G4double neutronMass = Neutron.getMass();
00186 G4double electronMass = 0.000511;
00187 G4double exnu = 0.;
00188 excitationEnergyGPN = 0.;
00189 excitationEnergyDTA = 0.;
00190
00191 if (atomicWeight > (neutronMass + 2.*electronMass))
00192 {
00193 G4int magic = ((G4int)(atomicNumber+0.1) == 82) ? 1 : 0;
00194 G4double ekin = Amin(Amax(incidentKineticEnergy, 0.1), 4.);
00195 G4double cfa = Amax(0.35 +((0.35 - 0.05)/2.3)*std::log(ekin), 0.15);
00196 exnu = 7.716*cfa*std::exp(-cfa);
00197 G4double atno = Amin(atomicWeight, 120.);
00198 cfa = ((atno - 1.)/120.) * std::exp(-(atno-1.)/120.);
00199 exnu = exnu * cfa;
00200 G4double fpdiv = Amax(1.-0.25*ekin*ekin, 0.5);
00201 G4double gfa = 2.*((atomicWeight-1.)/70.)
00202 * std::exp(-(atomicWeight-1.)/70.);
00203
00204 excitationEnergyGPN = exnu * fpdiv;
00205 excitationEnergyDTA = exnu - excitationEnergyGPN;
00206
00207 G4double ran1 = 0., ran2 = 0.;
00208 if (!magic)
00209 { ran1 = normal();
00210 ran2 = normal();
00211 }
00212 excitationEnergyGPN = Amax(excitationEnergyGPN*(1.+ran1*gfa),0.);
00213 excitationEnergyDTA = Amax(excitationEnergyDTA*(1.+ran2*gfa),0.);
00214 exnu = excitationEnergyGPN + excitationEnergyDTA;
00215 if(verboseLevel > 1) {
00216 G4cout << " NuclearExcitation: " << magic << " " << ekin
00217 << " " << cfa << " " << atno << " " << fpdiv << " "
00218 << gfa << " " << excitationEnergyGPN
00219 << " " << excitationEnergyDTA << G4endl;
00220 }
00221
00222 while (exnu >= incidentKineticEnergy)
00223 {
00224 excitationEnergyGPN *= (1. - 0.5*normal());
00225 excitationEnergyDTA *= (1. - 0.5*normal());
00226 exnu = excitationEnergyGPN + excitationEnergyDTA;
00227 }
00228 }
00229 return exnu;
00230 }
00231
00232 G4double
00233 G4HEInelastic::pmltpc(G4int npos, G4int nneg, G4int nzero, G4int n,
00234 G4double b, G4double c)
00235 {
00236 G4double expxu = 82.;
00237 G4double expxl = -expxu;
00238 G4int i;
00239 G4double npf = 0.0, nmf = 0.0, nzf = 0.0;
00240 for(i=2;i<=npos;i++) npf += std::log((G4double)i);
00241 for(i=2;i<=nneg;i++) nmf += std::log((G4double)i);
00242 for(i=2;i<=nzero;i++) nzf += std::log((G4double)i);
00243 G4double r = Amin(expxu,Amax(expxl,
00244 -(npos-nneg+nzero+b)*(npos-nneg+nzero+b)/(2*c*c*n*n)-npf-nmf-nzf));
00245 return std::exp(r);
00246 }
00247
00248
00249 G4int G4HEInelastic::Factorial(G4int n)
00250 {
00251 G4int result = 1;
00252 if (n < 0) G4Exception("G4HEInelastic::Factorial()", "HEP000",
00253 FatalException, "Negative factorial argument");
00254 while (n > 1) result *= n--;
00255 return result;
00256 }
00257
00258
00259 G4double G4HEInelastic::normal()
00260 {
00261 G4double ran = -6.0;
00262 for(G4int i=0; i<12; i++) ran += G4UniformRand();
00263 return ran;
00264 }
00265
00266
00267 G4int G4HEInelastic::Poisson(G4double x)
00268 {
00269 G4int i, iran = 0;
00270 G4double ran;
00271 if (x > 9.9) {
00272 iran = G4int( Amax( 0.0, x + normal() * std::sqrt( x ) ) );
00273 } else {
00274 G4int fivex = G4int(5.0 * x);
00275 if (fivex <= 0) {
00276 G4double p1 = x * std::exp( -x );
00277 G4double p2 = x * p1/2.;
00278 G4double p3 = x * p2/3.;
00279 ran = G4UniformRand();
00280 if (ran < p3) iran = 3;
00281 else if ( ran < p2 ) iran = 2;
00282 else if ( ran < p1 ) iran = 1;
00283 } else {
00284 G4double r = std::exp(-x);
00285 ran = G4UniformRand();
00286 if (ran > r) {
00287 G4double rrr;
00288 G4double rr = r;
00289 for (i = 1; i <= fivex; i++) {
00290 iran++;
00291 if (i > 5) rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((G4double)i)+i-0.9189385);
00292 else rrr = std::pow(x,i)*Factorial(i);
00293 rr += r * rrr;
00294 if (ran <= rr) break;
00295 }
00296 }
00297 }
00298 }
00299 return iran;
00300 }
00301
00302
00303 G4double
00304 G4HEInelastic::GammaRand( G4double avalue )
00305 {
00306 G4double ga = avalue -1.;
00307 G4double la = std::sqrt(2.*avalue - 1.);
00308 G4double ep = 1.570796327 + std::atan(ga/la);
00309 G4double ro = 1.570796327 - ep;
00310 G4double y = 1.;
00311 G4double xtrial;
00312 repeat:
00313 xtrial = ga + la * std::tan(ep*G4UniformRand() + ro);
00314 if(xtrial == 0.) goto repeat;
00315 y = std::log(1.+sqr((xtrial-ga)/la))+ga*std::log(xtrial/ga)-xtrial+ga;
00316 if(std::log(G4UniformRand()) > y) goto repeat;
00317 return xtrial;
00318 }
00319 G4double
00320 G4HEInelastic::Erlang( G4int mvalue )
00321 {
00322 G4double ran = G4UniformRand();
00323 G4double xtrial = 0.62666*std::log((1.+ran)/(1.-ran));
00324 if(G4UniformRand()<0.5) xtrial = -xtrial;
00325 return mvalue+xtrial*std::sqrt(G4double(mvalue));
00326 }
00327
00328 void
00329 G4HEInelastic::StrangeParticlePairProduction(
00330 const G4double availableEnergy,
00331 const G4double centerOfMassEnergy,
00332 G4HEVector pv[],
00333 G4int &vecLen,
00334 const G4HEVector& incidentParticle,
00335 const G4HEVector& targetParticle)
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 {
00346 static G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
00347 static G4double avkkb[] = {0.0015,0.0050,0.0120,0.0285,0.0525,0.0750,0.0975,
00348 0.1230,0.2800,0.3980,0.4950,0.5730};
00349 static G4double kkb[] = {0.2500,0.3750,0.5000,0.5625,0.6250,0.6875,0.7500,
00350 0.8750,1.0000};
00351 static G4double ky[] = {0.2000,0.3000,0.4000,0.5500,0.6250,0.7000,0.8000,
00352 0.8500,0.9000,0.9500,0.9750,1.0000};
00353 static G4int ipakkb[] = {10,13,10,11,10,12,11,11,11,12,12,11,12,12,
00354 11,13,12,13};
00355 static G4int ipaky[] = {18,10,18,11,18,12,20,10,20,11,20,12,21,10,
00356 21,11,21,12,22,10,22,11,22,12};
00357 static G4int ipakyb[] = {19,13,19,12,19,11,23,13,23,12,23,11,24,13,
00358 24,12,24,11,25,13,25,12,25,11};
00359 static G4double avky[] = {0.0050,0.0300,0.0640,0.0950,0.1150,0.1300,0.1450,
00360 0.1550,0.2000,0.2050,0.2100,0.2120};
00361 static G4double avnnb[] ={0.00001,0.0001,0.0006,0.0025,0.0100,0.0200,0.0400,
00362 0.0500,0.1200,0.1500,0.1800,0.2000};
00363
00364 G4int i, ibin, i3, i4;
00365 G4double avk, avy, avn, ran;
00366
00367 G4double protonMass = Proton.getMass();
00368 G4double sigmaMinusMass = SigmaMinus.getMass();
00369 G4int antiprotonCode = AntiProton.getCode();
00370 G4int antineutronCode = AntiNeutron.getCode();
00371 G4int antilambdaCode = AntiLambda.getCode();
00372
00373 G4double incidentMass = incidentParticle.getMass();
00374 G4int incidentCode = incidentParticle.getCode();
00375
00376 G4double targetMass = targetParticle.getMass();
00377
00378
00379
00380 if (vecLen <= 2) return;
00381
00382
00383
00384 i = 1;
00385 while ( (i<12) && (centerOfMassEnergy > avrs[i]) )i++;
00386 if ( i == 12 ) ibin = 11;
00387 else ibin = i;
00388
00389
00390
00391
00392 if( vecLen == 3 ) {
00393 i3 = 2;
00394 i4 = 3;
00395 }
00396 else
00397 {
00398 i4 = i3 = 2 + G4int( (vecLen-2)*G4UniformRand() );
00399 while ( i3 == i4 ) i4 = 2 + G4int( (vecLen-2)*G4UniformRand() );
00400 }
00401
00402
00403
00404 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
00405 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
00406 avk = std::exp(avk);
00407
00408 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
00409 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
00410 avy = std::exp(avy);
00411
00412 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
00413 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
00414 avn = std::exp(avn);
00415
00416 if ( avk+avy+avn <= 0.0 ) return;
00417
00418 if ( incidentMass < protonMass ) avy /= 2.0;
00419 avy += avk+avn;
00420 avk += avn;
00421
00422 ran = G4UniformRand();
00423 if ( ran < avn )
00424 {
00425 if ( availableEnergy < 2.0) return;
00426 if ( vecLen == 3 )
00427 {
00428 if ( G4UniformRand() < 0.5 )
00429 {
00430 pv[i3] = Neutron;;
00431 pv[vecLen++] = AntiNeutron;
00432 }
00433 else
00434 {
00435 pv[i3] = Proton;
00436 pv[vecLen++] = AntiProton;
00437 }
00438 }
00439 else
00440 {
00441 if ( G4UniformRand() < 0.5 )
00442 {
00443 pv[i3] = Neutron;
00444 pv[i4] = AntiNeutron;
00445 }
00446 else
00447 {
00448 pv[i3] = Proton;
00449 pv[i4] = AntiProton;
00450 }
00451 }
00452 }
00453 else if ( ran < avk )
00454 {
00455 if ( availableEnergy < 1.0) return;
00456 G4double ran1 = G4UniformRand();
00457 i = 0;
00458 while( (i<9) && (ran1>kkb[i]) )i++;
00459 if ( i == 9 ) return;
00460
00461
00462
00463
00464 switch( ipakkb[i*2] )
00465 {
00466 case 10: pv[i3] = KaonPlus; break;
00467 case 11: pv[i3] = KaonZeroShort;break;
00468 case 12: pv[i3] = KaonZeroLong; break;
00469 case 13: pv[i3] = KaonMinus; break;
00470 }
00471
00472 if( vecLen == 2 )
00473 {
00474 switch( ipakkb[i*2+1] )
00475 {
00476 case 10: pv[vecLen++] = KaonPlus; break;
00477 case 11: pv[vecLen++] = KaonZeroShort;break;
00478 case 12: pv[vecLen++] = KaonZeroLong; break;
00479 case 13: pv[vecLen++] = KaonMinus; break;
00480 }
00481 }
00482 else
00483 {
00484 switch( ipakkb[i*2+1] )
00485 {
00486 case 10: pv[i4] = KaonPlus; break;
00487 case 11: pv[i4] = KaonZeroShort;break;
00488 case 12: pv[i4] = KaonZeroLong; break;
00489 case 13: pv[i4] = KaonMinus; break;
00490 }
00491 }
00492 }
00493 else if ( ran < avy )
00494 {
00495 if( availableEnergy < 1.6) return;
00496 G4double ran1 = G4UniformRand();
00497 i = 0;
00498 while( (i<12) && (ran1>ky[i]) )i++;
00499 if ( i == 12 ) return;
00500 if ( (incidentMass<protonMass) || (G4UniformRand()<0.5) )
00501 {
00502
00503
00504
00505
00506
00507
00508
00509 switch( ipaky[i*2] )
00510 {
00511 case 18: pv[1] = Lambda; break;
00512 case 20: pv[1] = SigmaPlus; break;
00513 case 21: pv[1] = SigmaZero; break;
00514 case 22: pv[1] = SigmaMinus;break;
00515 }
00516 switch( ipaky[i*2+1] )
00517 {
00518 case 10: pv[i3] = KaonPlus; break;
00519 case 11: pv[i3] = KaonZeroShort;break;
00520 case 12: pv[i3] = KaonZeroLong; break;
00521 }
00522 }
00523 else
00524 {
00525
00526
00527
00528
00529
00530
00531 if( (incidentCode==antiprotonCode) || (incidentCode==antineutronCode) ||
00532 (incidentCode==antilambdaCode) || (incidentMass>sigmaMinusMass) )
00533 {
00534 switch( ipakyb[i*2] )
00535 {
00536 case 19:pv[0] = AntiLambda; break;
00537 case 23:pv[0] = AntiSigmaPlus; break;
00538 case 24:pv[0] = AntiSigmaZero; break;
00539 case 25:pv[0] = AntiSigmaMinus;break;
00540 }
00541 switch( ipakyb[i*2+1] )
00542 {
00543 case 11:pv[i3] = KaonZeroShort;break;
00544 case 12:pv[i3] = KaonZeroLong; break;
00545 case 13:pv[i3] = KaonMinus; break;
00546 }
00547 }
00548 else
00549 {
00550 switch( ipaky[i*2] )
00551 {
00552 case 18:pv[0] = Lambda; break;
00553 case 20:pv[0] = SigmaPlus; break;
00554 case 21:pv[0] = SigmaZero; break;
00555 case 22:pv[0] = SigmaMinus;break;
00556 }
00557 switch( ipaky[i*2+1] )
00558 {
00559 case 10: pv[i3] = KaonPlus; break;
00560 case 11: pv[i3] = KaonZeroShort;break;
00561 case 12: pv[i3] = KaonZeroLong; break;
00562 }
00563 }
00564 }
00565 }
00566 else
00567 return;
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577 incidentMass = incidentParticle.getMass();
00578 targetMass = targetParticle.getMass();
00579
00580 G4double energyCheck = centerOfMassEnergy-(incidentMass+targetMass);
00581 if (verboseLevel > 1) G4cout << "Particles produced: " ;
00582
00583 for ( i=0; i < vecLen; i++ )
00584 {
00585 energyCheck -= pv[i].getMass();
00586 if (verboseLevel > 1) G4cout << pv[i].getCode() << " " ;
00587 if( energyCheck < 0.0 )
00588 {
00589 if( i > 0 ) vecLen = --i;
00590 return;
00591 }
00592 }
00593 if (verboseLevel > 1) G4cout << G4endl;
00594 return;
00595 }
00596
00597 void
00598 G4HEInelastic::HighEnergyCascading(G4bool& successful,
00599 G4HEVector pv[],
00600 G4int& vecLen,
00601 G4double& excitationEnergyGNP,
00602 G4double& excitationEnergyDTA,
00603 const G4HEVector& incidentParticle,
00604 const G4HEVector& targetParticle,
00605 G4double atomicWeight,
00606 G4double atomicNumber)
00607 {
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 G4int protonCode = Proton.getCode();
00621 G4double protonMass = Proton.getMass();
00622 G4int neutronCode = Neutron.getCode();
00623 G4double neutronMass = Neutron.getMass();
00624 G4double kaonPlusMass = KaonPlus.getMass();
00625 G4int kaonPlusCode = KaonPlus.getCode();
00626 G4int kaonMinusCode = KaonMinus.getCode();
00627 G4int kaonZeroSCode = KaonZeroShort.getCode();
00628 G4int kaonZeroLCode = KaonZeroLong.getCode();
00629 G4int kaonZeroCode = KaonZero.getCode();
00630 G4int antiKaonZeroCode = AntiKaonZero.getCode();
00631 G4int pionPlusCode = PionPlus.getCode();
00632 G4int pionZeroCode = PionZero.getCode();
00633 G4int pionMinusCode = PionMinus.getCode();
00634 G4String mesonType = PionPlus.getType();
00635 G4String baryonType = Proton.getType();
00636 G4String antiBaryonType = AntiProton.getType();
00637
00638 G4double targetMass = targetParticle.getMass();
00639
00640 G4int incidentCode = incidentParticle.getCode();
00641 G4double incidentMass = incidentParticle.getMass();
00642 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00643 G4double incidentEnergy = incidentParticle.getEnergy();
00644 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
00645 G4String incidentType = incidentParticle.getType();
00646
00647 G4double incidentTOF = 0.;
00648
00649
00650
00651 G4int i, j, l;
00652
00653 if (verboseLevel > 1) G4cout << " G4HEInelastic::HighEnergyCascading "
00654 << G4endl;
00655 successful = false;
00656 if (incidentTotalMomentum < 25. + G4UniformRand()*25.) return;
00657
00658
00659
00660 G4bool annihilation = false;
00661 if (incidentCode < 0 && incidentType == antiBaryonType &&
00662 pv[0].getType() != antiBaryonType &&
00663 pv[1].getType() != antiBaryonType) {
00664 annihilation = true;
00665 }
00666
00667 G4double twsup[] = { 1., 1., 0.7, 0.5, 0.3, 0.2, 0.1, 0.0 };
00668
00669 if (annihilation) goto start;
00670 if (vecLen >= 8) goto start;
00671 if( incidentKineticEnergy < 1.) return;
00672 if( ( incidentCode == kaonPlusCode || incidentCode == kaonMinusCode
00673 || incidentCode == kaonZeroCode || incidentCode == antiKaonZeroCode
00674 || incidentCode == kaonZeroSCode || incidentCode == kaonZeroLCode )
00675 && ( G4UniformRand() < 0.5) ) goto start;
00676 if( G4UniformRand() > twsup[vecLen-1]) goto start;
00677 if( incidentKineticEnergy > (G4UniformRand()*200 + 50.) ) goto start;
00678 return;
00679
00680 start:
00681
00682 if (annihilation) {
00683
00684 G4double ekcor = Amax( 1., 1./incidentKineticEnergy);
00685 incidentKineticEnergy = 2*targetMass + incidentKineticEnergy*(1.+ekcor/atomicWeight);
00686 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00687 atomicWeight,
00688 atomicNumber,
00689 excitationEnergyGNP,
00690 excitationEnergyDTA);
00691 incidentKineticEnergy -= excitation;
00692 if (incidentKineticEnergy < excitationEnergyDTA) incidentKineticEnergy = 0.;
00693 incidentEnergy = incidentKineticEnergy + incidentMass;
00694 incidentTotalMomentum =
00695 std::sqrt( Amax(0., incidentEnergy*incidentEnergy - incidentMass*incidentMass));
00696 }
00697
00698 G4HEVector pTemp;
00699 for (i = 2; i < vecLen; i++) {
00700 j = Imin( vecLen-1, (G4int)(2. + G4UniformRand()*(vecLen - 2)));
00701 pTemp = pv[j];
00702 pv[j] = pv[i];
00703 pv[i] = pTemp;
00704 }
00705
00706
00707
00708
00709 if ((incidentCode==kaonPlusCode || incidentCode==kaonMinusCode ||
00710 incidentCode==kaonZeroCode || incidentCode==antiKaonZeroCode ||
00711 incidentCode==kaonZeroSCode || incidentCode==kaonZeroLCode)
00712 && (G4UniformRand() > 0.9) ) {
00713 pTemp = pv[1];
00714 pv[1] = pv[0];
00715 pv[0] = pTemp;
00716 }
00717
00718
00719
00720
00721
00722
00723
00724 G4int lead = 0;
00725 G4HEVector leadParticle;
00726 if ((incidentMass >= kaonPlusMass-0.05) &&
00727 (incidentCode != protonCode) && (incidentCode != neutronCode) ) {
00728 G4double pMass = pv[0].getMass();
00729 G4int pCode = pv[0].getCode();
00730 if ((pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
00731 && (pCode != neutronCode) ) {
00732 lead = pCode;
00733 leadParticle = pv[0];
00734 } else {
00735 pMass = pv[1].getMass();
00736 pCode = pv[1].getCode();
00737 if ((pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
00738 && (pCode != neutronCode) ) {
00739 lead = pCode;
00740 leadParticle = pv[1];
00741 }
00742 }
00743 }
00744
00745
00746
00747
00748 G4HEVector pvI = incidentParticle;
00749 pvI.setSide( 1 );
00750
00751 G4HEVector pvT = targetParticle;
00752 pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
00753 pvT.setSide( -1 );
00754 pvT.setTOF( -1.);
00755
00756 G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass)+sqr(targetMass)
00757 +2.0*targetMass*incidentEnergy );
00758
00759 G4double tavai1 = centerOfMassEnergy/2.0 - incidentMass;
00760 G4double tavai2 = centerOfMassEnergy/2.0 - targetMass;
00761
00762
00763
00764
00765 G4int ntb = 1;
00766 for (i = 0; i < vecLen; i++) {
00767 if (i == 0) pv[i].setSide(1);
00768 else if (i == 1) pv[i].setSide(-1);
00769 else {
00770 if (G4UniformRand() < 0.5) {
00771 pv[i].setSide(-1);
00772 ntb++;
00773 } else {
00774 pv[i].setSide(1);
00775 }
00776 }
00777 pv[i].setTOF(incidentTOF);
00778 }
00779
00780 G4double tb = 2. * ntb;
00781 if (centerOfMassEnergy < (2. + G4UniformRand()))
00782 tb = (2. * ntb + vecLen)/2.;
00783
00784 if (verboseLevel > 1) {
00785 G4cout << " pv Vector after Randomization " << vecLen << G4endl;
00786 pvI.Print(-1);
00787 pvT.Print(-1);
00788 for (i = 0; i < vecLen; i++) pv[i].Print(i);
00789 }
00790
00791
00792
00793
00794
00795
00796 G4double ss, xtarg, ran;
00797 ss = centerOfMassEnergy*centerOfMassEnergy;
00798 G4double afc;
00799 afc = Amin(0.5, 0.312 + 0.200 * std::log(std::log(ss))+ std::pow(ss,1.5)/6000.0);
00800 xtarg = Amax(0.01, afc * (std::pow(atomicWeight, 0.33) - 1.0) * tb);
00801 G4int nstran = Poisson( 0.03*xtarg);
00802 G4int momentumBin = 0;
00803 G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.5, 0.5 };
00804 G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
00805 while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin])) momentumBin++;
00806 momentumBin = Imin(5, momentumBin);
00807 G4double xpnhmf = Amax(0.01,xtarg*nucsup[momentumBin]);
00808 G4double xshhmf = Amax(0.01,xtarg - xpnhmf);
00809 G4double rshhmf = 0.25*xshhmf;
00810 G4double rpnhmf = 0.81*xpnhmf;
00811 G4double xhmf=0;
00812 if (verboseLevel > 1)
00813 G4cout << "xtarg= " << xtarg << " xpnhmf = " << xpnhmf << G4endl;
00814
00815 G4int nshhmf, npnhmf;
00816 if (rshhmf > 1.1) {
00817 rshhmf = xshhmf/(rshhmf-1.);
00818 if (rshhmf <= 20.)
00819 xhmf = GammaRand( rshhmf );
00820 else
00821 xhmf = Erlang( G4int(rshhmf+0.5) );
00822 xshhmf *= xhmf/rshhmf;
00823 }
00824 nshhmf = Poisson( xshhmf );
00825 if(verboseLevel > 1)
00826 G4cout << "xshhmf = " << xshhmf << " xhmf = " << xhmf
00827 << " rshhmf = " << rshhmf << G4endl;
00828
00829 if (rpnhmf > 1.1)
00830 {
00831 rpnhmf = xpnhmf/(rpnhmf -1.);
00832 if (rpnhmf <= 20.)
00833 xhmf = GammaRand( rpnhmf );
00834 else
00835 xhmf = Erlang( G4int(rpnhmf+0.5) );
00836 xpnhmf *= xhmf/rpnhmf;
00837 }
00838 npnhmf = Poisson( xpnhmf );
00839 if(verboseLevel > 1)
00840 G4cout << "nshhmf = " << nshhmf << " npnhmf = " << npnhmf
00841 << " nstran = " << nstran << G4endl;
00842
00843 G4int ntarg = nshhmf + npnhmf + nstran;
00844
00845 G4int targ = 0;
00846
00847 while (npnhmf > 0)
00848 {
00849 if ( G4UniformRand() > (1. - atomicNumber/atomicWeight))
00850 pv[vecLen] = Proton;
00851 else
00852 pv[vecLen] = Neutron;
00853 targ++;
00854 pv[vecLen].setSide( -2 );
00855 pv[vecLen].setFlag( true );
00856 pv[vecLen].setTOF( incidentTOF );
00857 vecLen++;
00858 npnhmf--;
00859 }
00860 while (nstran > 0)
00861 {
00862 ran = G4UniformRand();
00863 if (ran < 0.14) pv[vecLen] = Lambda;
00864 else if (ran < 0.20) pv[vecLen] = SigmaZero;
00865 else if (ran < 0.43) pv[vecLen] = KaonPlus;
00866 else if (ran < 0.66) pv[vecLen] = KaonZero;
00867 else if (ran < 0.89) pv[vecLen] = AntiKaonZero;
00868 else pv[vecLen] = KaonMinus;
00869 if (G4UniformRand() > 0.2)
00870 {
00871 pv[vecLen].setSide( -2 );
00872 pv[vecLen].setFlag( true );
00873 }
00874 else
00875 {
00876 pv[vecLen].setSide( 1 );
00877 pv[vecLen].setFlag( false );
00878 ntarg--;
00879 }
00880 pv[vecLen].setTOF( incidentTOF );
00881 vecLen++;
00882 nstran--;
00883 }
00884 while (nshhmf > 0)
00885 {
00886 ran = G4UniformRand();
00887 if( ran < 0.33333 )
00888 pv[vecLen] = PionPlus;
00889 else if( ran < 0.66667 )
00890 pv[vecLen] = PionZero;
00891 else
00892 pv[vecLen] = PionMinus;
00893 if (G4UniformRand() > 0.2)
00894 {
00895 pv[vecLen].setSide( -2 );
00896 pv[vecLen].setFlag( true );
00897 }
00898 else
00899 {
00900 pv[vecLen].setSide( 1 );
00901 pv[vecLen].setFlag( false );
00902 ntarg--;
00903 }
00904 pv[vecLen].setTOF( incidentTOF );
00905 vecLen++;
00906 nshhmf--;
00907 }
00908
00909
00910
00911
00912 G4int is, iskip, iavai1;
00913 if (vecLen <= 1) return;
00914
00915 tavai1 = centerOfMassEnergy/2.;
00916 iavai1 = 0;
00917
00918 for (i = 0; i < vecLen; i++)
00919 {
00920 if (pv[i].getSide() > 0)
00921 {
00922 tavai1 -= pv[i].getMass();
00923 iavai1++;
00924 }
00925 }
00926 if ( iavai1 == 0) return;
00927
00928 while (tavai1 <= 0.0) {
00929
00930 iskip = G4int(G4UniformRand()*iavai1) + 1;
00931 is = 0;
00932 for (i = vecLen-1; i >= 0; i--) {
00933 if (pv[i].getSide() > 0) {
00934 if (++is == iskip) {
00935 tavai1 += pv[i].getMass();
00936 iavai1--;
00937 if (i != vecLen-1) {
00938 for (j = i; j < vecLen; j++) pv[j] = pv[j+1];
00939 }
00940 if (--vecLen == 0) return;
00941 break;
00942 }
00943 }
00944 }
00945 }
00946
00947 tavai2 = (targ+1)*centerOfMassEnergy/2.;
00948 G4int iavai2 = 0;
00949
00950 for (i = 0; i < vecLen; i++)
00951 {
00952 if (pv[i].getSide() < 0)
00953 {
00954 tavai2 -= pv[i].getMass();
00955 iavai2++;
00956 }
00957 }
00958 if (iavai2 == 0) return;
00959
00960 while( tavai2 <= 0.0 )
00961 {
00962 iskip = G4int(G4UniformRand()*iavai2) + 1;
00963 is = 0;
00964 for( i = vecLen-1; i >= 0; i-- )
00965 {
00966 if( pv[i].getSide() < 0 )
00967 {
00968 if( ++is == iskip )
00969 {
00970 tavai2 += pv[i].getMass();
00971 iavai2--;
00972 if (pv[i].getSide() == -2) ntarg--;
00973 if (i != vecLen-1)
00974 {
00975 for( j=i; j<vecLen; j++)
00976 {
00977 pv[j] = pv[j+1];
00978 }
00979 }
00980 if (--vecLen == 0) return;
00981 break;
00982 }
00983 }
00984 }
00985 }
00986
00987 if (verboseLevel > 1) {
00988 G4cout << " pv Vector after Energy checks "
00989 << vecLen << " " << tavai1 << " " << iavai1 << " " << tavai2
00990 << " " << iavai2 << " " << ntarg << G4endl;
00991 pvI.Print(-1);
00992 pvT.Print(-1);
00993 for (i=0; i < vecLen ; i++) pv[i].Print(i);
00994 }
00995
00996
00997
00998 G4HEVector* pvmx = new G4HEVector [10];
00999
01000 pvmx[0].setMass( incidentMass );
01001 pvmx[0].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
01002 pvmx[1].setMass( protonMass);
01003 pvmx[1].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
01004 pvmx[3].setMass( protonMass*(1+targ));
01005 pvmx[3].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
01006 pvmx[4].setZero();
01007 pvmx[5].setZero();
01008 pvmx[7].setZero();
01009 pvmx[8].setZero();
01010 pvmx[8].setMomentum( 1.0, 0.0 );
01011 pvmx[2].Add( pvmx[0], pvmx[1] );
01012 pvmx[3].Add( pvmx[3], pvmx[0] );
01013 pvmx[0].Lor( pvmx[0], pvmx[2] );
01014 pvmx[1].Lor( pvmx[1], pvmx[2] );
01015
01016 if (verboseLevel > 1) {
01017 G4cout << " General Vectors after Definition " << G4endl;
01018 for (i=0; i<10; i++) pvmx[i].Print(i);
01019 }
01020
01021
01022
01023
01024
01025 G4double dndl[20];
01026 G4double binl[20];
01027 G4double pvMass(0), pvEnergy(0);
01028 G4int pvCode;
01029 G4double aspar, pt, phi, et, xval;
01030 G4double ekin = 0.;
01031 G4double ekin1 = 0.;
01032 G4double ekin2 = 0.;
01033 G4int npg = 0;
01034 G4double rmg0 = 0.;
01035 G4int targ1 = 0;
01036 phi = G4UniformRand()*twopi;
01037
01038 for (i = vecLen-1; i >= 0; i--) {
01039
01040 if( i == 1)
01041 {
01042 if ( (pv[i].getMass() > neutronMass + 0.05) && (G4UniformRand() < 0.2))
01043 {
01044 if(++npg < 19)
01045 {
01046 pv[i].setSide(-3);
01047 rmg0 += pv[i].getMass();
01048 targ++;
01049 continue;
01050 }
01051 }
01052 else if ( pv[i].getMass() > protonMass - 0.05)
01053 {
01054 if(++npg < 19)
01055 {
01056 pv[i].setSide(-3);
01057 rmg0 += pv[i].getMass();
01058 targ++;
01059 continue;
01060 }
01061 }
01062 }
01063 if( pv[i].getSide() == -2)
01064 {
01065 if ( pv[i].getName() == "Proton" || pv[i].getName() == "Neutron")
01066 {
01067 if( ++npg < 19 )
01068 {
01069 pv[i].setSide( -3 );
01070 rmg0 += pv[i].getMass();
01071 targ1++;
01072 continue;
01073 }
01074 }
01075 }
01076
01077
01078
01079
01080 G4double maspar[] = { 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.20, 0.10};
01081 G4double bp[] = { 4.00, 2.50, 2.20, 3.00, 3.00, 1.70, 3.50, 3.50};
01082 G4double ptex[] = { 1.70, 1.70, 1.50, 1.70, 1.40, 1.20, 1.70, 1.20};
01083
01084
01085
01086
01087
01088 pvMass = pv[i].getMass();
01089 j = 2;
01090 if (pv[i].getType() == mesonType ) j = 1;
01091 if ( pv[i].getMass() < 0.4 ) j = 0;
01092 if ( i <= 1 ) j += 3;
01093 if (pv[i].getSide() <= -2) j = 6;
01094 if (j == 6 && (pv[i].getType() == baryonType || pv[i].getType() == antiBaryonType)) j = 7;
01095 pt = std::sqrt(std::pow(-std::log(1.-G4UniformRand())/bp[j],ptex[j]));
01096 if(pt<0.05) pt = Amax(0.001, 0.3*G4UniformRand());
01097 aspar = maspar[j];
01098 phi = G4UniformRand()*twopi;
01099 pv[i].setMomentum( pt*std::cos(phi), pt*std::sin(phi) );
01100
01101 for( j=0; j<20; j++ ) binl[j] = j/(19.*pt);
01102
01103 if( pv[i].getSide() > 0 )
01104 et = pvmx[0].getEnergy();
01105 else
01106 et = pvmx[1].getEnergy();
01107
01108 dndl[0] = 0.0;
01109
01110
01111
01112 G4int outerCounter = 0, innerCounter = 0;
01113 G4bool eliminateThisParticle = true;
01114 G4bool resetEnergies = true;
01115 while( ++outerCounter < 3 )
01116 {
01117 for( l=1; l<20; l++ )
01118 {
01119 xval = (binl[l]+binl[l-1])/2.;
01120 if( xval > 1./pt )
01121 dndl[l] = dndl[l-1];
01122 else
01123 dndl[l] = dndl[l-1] +
01124 aspar/std::sqrt( std::pow((1.+aspar*xval*aspar*xval),3) ) *
01125 (binl[l]-binl[l-1]) * et /
01126 std::sqrt( pt*xval*et*pt*xval*et + pt*pt + pvMass*pvMass );
01127 }
01128
01129
01130
01131 innerCounter = 0;
01132 while( ++innerCounter < 7 )
01133 {
01134 l = 1;
01135 ran = G4UniformRand()*dndl[19];
01136 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
01137 l = Imin( 19, l );
01138 xval = Amin( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
01139 if( pv[i].getSide() < 0 ) xval *= -1.;
01140 pv[i].setMomentumAndUpdate( xval*et );
01141 pvEnergy = pv[i].getEnergy();
01142 if( pv[i].getSide() > 0 )
01143 {
01144 if ( i < 2 )
01145 {
01146 ekin = tavai1 - ekin1;
01147 if (ekin < 0.) ekin = 0.04*std::fabs(normal());
01148 G4double pp1 = pv[i].Length();
01149 if (pp1 >= 1.e-6)
01150 {
01151 G4double pp = std::sqrt(ekin*(ekin+2*pvMass));
01152 pp = Amax(0., pp*pp - pt*pt);
01153 pp = std::sqrt(pp)*pv[i].getSide()/std::fabs(G4double(pv[i].getSide()));
01154 pv[i].setMomentumAndUpdate( pp );
01155 }
01156 else
01157 {
01158 pv[i].setMomentum(0.,0.,0.);
01159 pv[i].setKineticEnergyAndUpdate( ekin);
01160 }
01161 pvmx[4].Add( pvmx[4], pv[i]);
01162 outerCounter = 2;
01163 resetEnergies = false;
01164 eliminateThisParticle = false;
01165 break;
01166 }
01167 else if( (ekin1+pvEnergy-pvMass) < 0.95*tavai1 )
01168 {
01169 pvmx[4].Add( pvmx[4], pv[i] );
01170 ekin1 += pvEnergy - pvMass;
01171 pvmx[6].Add( pvmx[4], pvmx[5] );
01172 pvmx[6].setMomentum( 0.0 );
01173 outerCounter = 2;
01174 eliminateThisParticle = false;
01175 resetEnergies = false;
01176 break;
01177 }
01178 if( innerCounter > 5 ) break;
01179
01180 if( tavai2 >= pvMass )
01181 {
01182 pv[i].setSide( -1 );
01183 tavai1 += pvMass;
01184 tavai2 -= pvMass;
01185 iavai2++;
01186 }
01187 }
01188 else
01189 {
01190 xval = Amin(0.999,0.95+0.05*targ/20.0);
01191 if( (ekin2+pvEnergy-pvMass) < xval*tavai2 )
01192 {
01193 pvmx[5].Add( pvmx[5], pv[i] );
01194 ekin2 += pvEnergy - pvMass;
01195 pvmx[6].Add( pvmx[4], pvmx[5] );
01196 pvmx[6].setMomentum( 0.0 );
01197 outerCounter = 2;
01198 eliminateThisParticle = false;
01199 resetEnergies = false;
01200 break;
01201 }
01202 if( innerCounter > 5 )break;
01203
01204 if( tavai1 >= pvMass )
01205 {
01206 pv[i].setSide( 1 );
01207 tavai1 -= pvMass;
01208 tavai2 += pvMass;
01209 iavai2--;
01210 }
01211 }
01212 pv[i].setMomentum( pv[i].getMomentum().x() * 0.9,
01213 pv[i].getMomentum().y() * 0.9);
01214 pt *= 0.9;
01215 dndl[19] *= 0.9;
01216 }
01217
01218 if (resetEnergies)
01219 {
01220 if (verboseLevel > 1) {
01221 G4cout << " Reset energies for index " << i << " "
01222 << ekin1 << " " << tavai1 << G4endl;
01223 pv[i].Print(i);
01224 }
01225 ekin1 = 0.0;
01226 ekin2 = 0.0;
01227 pvmx[4].setZero();
01228 pvmx[5].setZero();
01229
01230 for( l=i+1; l < vecLen; l++ )
01231 {
01232 if( (pv[l].getMass() < protonMass) || (pv[l].getSide() > 0) )
01233 {
01234 pvEnergy = pv[l].getMass() + 0.95*pv[l].getKineticEnergy();
01235 pv[l].setEnergyAndUpdate( pvEnergy );
01236 if( pv[l].getSide() > 0)
01237 {
01238 ekin1 += pv[l].getKineticEnergy();
01239 pvmx[4].Add( pvmx[4], pv[l] );
01240 }
01241 else
01242 {
01243 ekin2 += pv[l].getKineticEnergy();
01244 pvmx[5].Add( pvmx[5], pv[l] );
01245 }
01246 }
01247 }
01248 }
01249 }
01250
01251 if (eliminateThisParticle) {
01252
01253 if (verboseLevel > 1) {
01254 G4cout << " Eliminate particle index " << i << G4endl;
01255 pv[i].Print(i);
01256 }
01257 if (i != vecLen-1) {
01258
01259 for (j = i; j < vecLen-1; j++) pv[j] = pv[j+1];
01260 }
01261 vecLen--;
01262
01263 if (vecLen < 2) {
01264 delete [] pvmx;
01265 return;
01266 }
01267 i++;
01268 pvmx[6].Add( pvmx[4], pvmx[5] );
01269 pvmx[6].setMomentum( 0.0 );
01270 }
01271 }
01272
01273 if (verboseLevel > 1) {
01274 G4cout << " pv Vector after lambda fragmentation " << vecLen << G4endl;
01275 pvI.Print(-1);
01276 pvT.Print(-1);
01277 for (i=0; i < vecLen ; i++) pv[i].Print(i);
01278 for (i=0; i < 10; i++) pvmx[i].Print(i);
01279 }
01280
01281
01282
01283 G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
01284 G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
01285
01286 if (npg > 0) {
01287 G4double rmg = rmg0;
01288 if (npg > 1) {
01289 G4int npg1 = npg-1;
01290 if (npg1 >4) npg1 = 4;
01291 rmg += std::pow( -std::log(1.-G4UniformRand()), cpar[npg1])/gpar[npg1];
01292 }
01293 G4double ga = 1.2;
01294 G4double ekit1 = 0.04, ekit2 = 0.6;
01295 if (incidentKineticEnergy < 5.) {
01296 ekit1 *= sqr(incidentKineticEnergy)/25.;
01297 ekit2 *= sqr(incidentKineticEnergy)/25.;
01298 }
01299 G4double avalue = (1.-ga)/(std::pow(ekit2,1.-ga)-std::pow(ekit1,1.-ga));
01300 for (i = 0; i < vecLen; i++) {
01301 if (pv[i].getSide() == -3) {
01302 G4double ekit = std::pow(G4UniformRand()*(1.-ga)/avalue + std::pow(ekit1,1.-ga), 1./(1.-ga) );
01303 G4double cost = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
01304 G4double sint = std::sqrt(1. - cost*cost);
01305 G4double pphi = twopi*G4UniformRand();
01306 G4double pp = std::sqrt(ekit*(ekit+2*pv[i].getMass()));
01307 pv[i].setMomentum(pp*sint*std::sin(pphi),
01308 pp*sint*std::cos(pphi),
01309 pp*cost);
01310 pv[i].Lor( pv[i], pvmx[2] );
01311 pvmx[5].Add( pvmx[5], pv[i] );
01312 }
01313 }
01314 }
01315
01316 if (vecLen <= 2) {
01317 successful = false;
01318 delete [] pvmx;
01319 return;
01320 }
01321
01322
01323
01324 targ = 0;
01325 for( i=0; i < vecLen; i++ )
01326 {
01327 if( pv[i].getType() == baryonType )targ++;
01328 if( pv[i].getType() == antiBaryonType )targ--;
01329 if(verboseLevel > 1) pv[i].Print(i);
01330 pv[i].Lor( pv[i], pvmx[1] );
01331 if(verboseLevel > 1) pv[i].Print(i);
01332 }
01333 if ( targ <1) targ = 1;
01334
01335 G4bool dum=0;
01336 if( lead )
01337 {
01338 for( i=0; i<vecLen; i++ )
01339 {
01340 if( pv[i].getCode() == lead )
01341 {
01342 dum = false;
01343 break;
01344 }
01345 }
01346 if( dum )
01347 {
01348 i = 0;
01349
01350 if( ( (leadParticle.getType() == baryonType ||
01351 leadParticle.getType() == antiBaryonType)
01352 && (pv[1].getType() == baryonType ||
01353 pv[1].getType() == antiBaryonType))
01354 || ( (leadParticle.getType() == mesonType)
01355 && (pv[1].getType() == mesonType)))
01356 {
01357 i = 1;
01358 }
01359 ekin = pv[i].getKineticEnergy();
01360 pv[i] = leadParticle;
01361 if( pv[i].getFlag() )
01362 pv[i].setTOF( -1.0 );
01363 else
01364 pv[i].setTOF( 1.0 );
01365 pv[i].setKineticEnergyAndUpdate( ekin );
01366 }
01367 }
01368
01369 pvmx[3].setMass( incidentMass);
01370 pvmx[3].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
01371
01372 G4double ekin0 = pvmx[3].getKineticEnergy();
01373
01374 pvmx[4].setMass( protonMass * targ);
01375 pvmx[4].setEnergy( protonMass * targ);
01376 pvmx[4].setKineticEnergy(0.);
01377 pvmx[4].setMomentum(0., 0., 0.);
01378 ekin = pvmx[3].getEnergy() + pvmx[4].getEnergy();
01379
01380 pvmx[5].Add( pvmx[3], pvmx[4] );
01381 pvmx[3].Lor( pvmx[3], pvmx[5] );
01382 pvmx[4].Lor( pvmx[4], pvmx[5] );
01383
01384 G4double tecm = pvmx[3].getEnergy() + pvmx[4].getEnergy();
01385
01386 pvmx[7].setZero();
01387
01388 ekin1 = 0.0;
01389 G4double teta, wgt;
01390
01391 for( i=0; i < vecLen; i++ )
01392 {
01393 pvmx[7].Add( pvmx[7], pv[i] );
01394 ekin1 += pv[i].getKineticEnergy();
01395 ekin -= pv[i].getMass();
01396 }
01397
01398 if( vecLen > 1 && vecLen < 19 )
01399 {
01400 G4bool constantCrossSection = true;
01401 G4HEVector pw[19];
01402 for(i=0; i<vecLen; i++) pw[i] = pv[i];
01403 wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
01404 ekin = 0.0;
01405 for( i=0; i < vecLen; i++ )
01406 {
01407 pvmx[6].setMass( pw[i].getMass());
01408 pvmx[6].setMomentum( pw[i].getMomentum() );
01409 pvmx[6].SmulAndUpdate( pvmx[6], 1. );
01410 pvmx[6].Lor( pvmx[6], pvmx[4] );
01411 ekin += pvmx[6].getKineticEnergy();
01412 }
01413 teta = pvmx[7].Ang( pvmx[3] );
01414 if (verboseLevel > 1)
01415 G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
01416 << " " << ekin1 << " " << ekin << G4endl;
01417 }
01418
01419 if( ekin1 != 0.0 )
01420 {
01421 pvmx[6].setZero();
01422 wgt = ekin/ekin1;
01423 ekin1 = 0.;
01424 for( i=0; i < vecLen; i++ )
01425 {
01426 pvMass = pv[i].getMass();
01427 ekin = pv[i].getKineticEnergy() * wgt;
01428 pv[i].setKineticEnergyAndUpdate( ekin );
01429 ekin1 += ekin;
01430 pvmx[6].Add( pvmx[6], pv[i] );
01431 }
01432 teta = pvmx[6].Ang( pvmx[3] );
01433 if (verboseLevel > 1) {
01434 G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
01435 << ekin1 << G4endl;
01436 incidentParticle.Print(0);
01437 targetParticle.Print(1);
01438 for(i=0;i<vecLen;i++) pv[i].Print(i);
01439 }
01440 }
01441
01442
01443
01444 G4double ry = G4UniformRand();
01445 G4double rz = G4UniformRand();
01446 G4double rx = twopi*rz;
01447 G4double a1 = std::sqrt(-2.0*std::log(ry));
01448 G4double rantarg1 = a1*std::cos(rx)*0.02*targ/G4double(vecLen);
01449 G4double rantarg2 = a1*std::sin(rx)*0.02*targ/G4double(vecLen);
01450
01451 for (i = 0; i < vecLen; i++)
01452 pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
01453 pv[i].getMomentum().y()+rantarg2 );
01454
01455 if (verboseLevel > 1) {
01456 pvmx[6].setZero();
01457 for (i = 0; i < vecLen; i++) pvmx[6].Add( pvmx[6], pv[i] );
01458 teta = pvmx[6].Ang( pvmx[3] );
01459 G4cout << " After smearing " << teta << G4endl;
01460 }
01461
01462
01463
01464
01465
01466
01467
01468
01469 G4double dekin = 0.0;
01470 G4int npions = 0;
01471 G4double ek1 = 0.0;
01472 G4double alekw, xxh;
01473 G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
01474 G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.00};
01475 G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
01476
01477 if (verboseLevel > 1)
01478 G4cout << " Rotation in Direction of primary particle (Defs1)" << G4endl;
01479
01480 for (i = 0; i < vecLen; i++) {
01481 if(verboseLevel > 1) pv[i].Print(i);
01482 pv[i].Defs1( pv[i], pvI );
01483 if(verboseLevel > 1) pv[i].Print(i);
01484 if (atomicWeight > 1.5) {
01485 ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa*( 1. + 0.5*normal()));
01486 alekw = std::log( incidentKineticEnergy );
01487 xxh = 1.;
01488 if (incidentCode == pionPlusCode || incidentCode == pionMinusCode) {
01489 if (pv[i].getCode() == pionZeroCode) {
01490 if (G4UniformRand() < std::log(atomicWeight)) {
01491 if (alekw > alem[0]) {
01492 G4int jmax = 1;
01493 for (j = 1; j < 8; j++) {
01494 if (alekw < alem[j]) {
01495 jmax = j;
01496 break;
01497 }
01498 }
01499 xxh = (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alekw
01500 + val0[jmax-1] - (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alem[jmax-1];
01501 xxh = 1. - xxh;
01502 }
01503 }
01504 }
01505 }
01506 dekin += ekin*(1.-xxh);
01507 ekin *= xxh;
01508 pv[i].setKineticEnergyAndUpdate( ekin );
01509 pvCode = pv[i].getCode();
01510 if ((pvCode == pionPlusCode) ||
01511 (pvCode == pionMinusCode) ||
01512 (pvCode == pionZeroCode)) {
01513 npions += 1;
01514 ek1 += ekin;
01515 }
01516 }
01517 }
01518
01519 if ( (ek1 > 0.0) && (npions > 0) ) {
01520 dekin = 1.+dekin/ek1;
01521 for (i = 0; i < vecLen; i++) {
01522 pvCode = pv[i].getCode();
01523 if ((pvCode == pionPlusCode) ||
01524 (pvCode == pionMinusCode) ||
01525 (pvCode == pionZeroCode)) {
01526 ekin = Amax(1.0e-6, pv[i].getKineticEnergy() * dekin);
01527 pv[i].setKineticEnergyAndUpdate( ekin );
01528 }
01529 }
01530 }
01531
01532 if (verboseLevel > 1) {
01533 G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
01534 incidentParticle.Print(0);
01535 targetParticle.Print(1);
01536 for (i = 0; i < vecLen; i++) pv[i].Print(i);
01537 }
01538
01539
01540
01541
01542
01543 if (verboseLevel > 1)
01544 G4cout << " Evaporation : " << atomicWeight << " "
01545 << excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
01546
01547 G4double sprob = 0.;
01548 if (incidentKineticEnergy > 5.)
01549
01550 sprob = Amin(1., 0.000314*atomicWeight*std::log(incidentKineticEnergy-4.));
01551 if( atomicWeight > 1.5 && G4UniformRand() > sprob )
01552 {
01553
01554 G4double cost, sint, pp, eka;
01555 G4int spall(0), nbl(0);
01556
01557
01558
01559 if( excitationEnergyGNP >= 0.001 )
01560 {
01561
01562
01563
01564 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyGNP/
01565 (excitationEnergyGNP+excitationEnergyDTA));
01566 if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
01567 if (verboseLevel > 1)
01568 G4cout << " evaporation " << targ << " " << nbl << " "
01569 << sprob << G4endl;
01570 spall = targ;
01571 if( nbl > 0)
01572 {
01573 ekin = (excitationEnergyGNP)/nbl;
01574 ekin2 = 0.0;
01575 for( i=0; i<nbl; i++ )
01576 {
01577 if( G4UniformRand() < sprob )
01578 {
01579 if(verboseLevel > 1) G4cout << " Particle skipped " << G4endl;
01580 continue;
01581 }
01582 if( ekin2 > excitationEnergyGNP) break;
01583 ran = G4UniformRand();
01584 ekin1 = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
01585 if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
01586 ekin2 += ekin1;
01587 if( ekin2 > excitationEnergyGNP)
01588 ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
01589 if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
01590 pv[vecLen] = Proton;
01591 else
01592 pv[vecLen] = Neutron;
01593 spall++;
01594 cost = G4UniformRand() * 2.0 - 1.0;
01595 sint = std::sqrt(std::fabs(1.0-cost*cost));
01596 phi = twopi * G4UniformRand();
01597 pv[vecLen].setFlag( true );
01598 pv[vecLen].setSide( -4 );
01599 pv[vecLen].setTOF( 1.0 );
01600 pvMass = pv[vecLen].getMass();
01601 pvEnergy = ekin1 + pvMass;
01602 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
01603 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
01604 pp*sint*std::cos(phi),
01605 pp*cost );
01606 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
01607 vecLen++;
01608 }
01609 if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
01610 {
01611 G4int ika, kk = 0;
01612 eka = incidentKineticEnergy;
01613 if( eka > 1.0 )eka *= eka;
01614 eka = Amax( 0.1, eka );
01615 ika = G4int(3.6*std::exp((atomicNumber*atomicNumber
01616 /atomicWeight-35.56)/6.45)/eka);
01617 if( ika > 0 )
01618 {
01619 for( i=(vecLen-1); i>=0; i-- )
01620 {
01621 if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
01622 {
01623 pTemp = pv[i];
01624 pv[i].setDefinition("Neutron");
01625 pv[i].setMomentumAndUpdate(pTemp.getMomentum());
01626 if (verboseLevel > 1) pv[i].Print(i);
01627 if( ++kk > ika ) break;
01628 }
01629 }
01630 }
01631 }
01632 }
01633 }
01634
01635
01636
01637
01638 if( excitationEnergyDTA >= 0.001 )
01639 {
01640 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyDTA
01641 /(excitationEnergyGNP+excitationEnergyDTA));
01642
01643
01644
01645 if (verboseLevel > 1)
01646 G4cout << " evaporation " << targ << " " << nbl << " "
01647 << sprob << G4endl;
01648 if( nbl > 0 )
01649 {
01650 ekin = excitationEnergyDTA/nbl;
01651 ekin2 = 0.0;
01652 for( i=0; i<nbl; i++ )
01653 {
01654 if( G4UniformRand() < sprob )
01655 {
01656 if(verboseLevel > 1) G4cout << " Particle skipped " << G4endl;
01657 continue;
01658 }
01659 if( ekin2 > excitationEnergyDTA) break;
01660 ran = G4UniformRand();
01661 ekin1 = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
01662 if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
01663 ekin2 += ekin1;
01664 if( ekin2 > excitationEnergyDTA)
01665 ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
01666 cost = G4UniformRand()*2.0 - 1.0;
01667 sint = std::sqrt(std::fabs(1.0-cost*cost));
01668 phi = twopi*G4UniformRand();
01669 ran = G4UniformRand();
01670 if( ran <= 0.60 )
01671 pv[vecLen] = Deuteron;
01672 else if (ran <= 0.90)
01673 pv[vecLen] = Triton;
01674 else
01675 pv[vecLen] = Alpha;
01676 spall += (int)(pv[vecLen].getMass() * 1.066);
01677 if( spall > atomicWeight ) break;
01678 pv[vecLen].setFlag( true );
01679 pv[vecLen].setSide( -4 );
01680 pvMass = pv[vecLen].getMass();
01681 pv[vecLen].setTOF( 1.0 );
01682 pvEnergy = pvMass + ekin1;
01683 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
01684 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
01685 pp*sint*std::cos(phi),
01686 pp*cost );
01687 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
01688 vecLen++;
01689 }
01690 }
01691 }
01692 }
01693 if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
01694 {
01695 for( i=0; i<vecLen; i++ )
01696 {
01697 G4double etb = pv[i].getKineticEnergy();
01698 if( etb >= incidentKineticEnergy )
01699 pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
01700 }
01701 }
01702
01703 if(verboseLevel > 1)
01704 { G4cout << "Call TuningOfHighEnergyCacsading vecLen = " << vecLen << G4endl;
01705 incidentParticle.Print(0);
01706 targetParticle.Print(1);
01707 for (i=0; i<vecLen; i++) pv[i].Print(i);
01708 }
01709
01710 TuningOfHighEnergyCascading( pv, vecLen,
01711 incidentParticle, targetParticle,
01712 atomicWeight, atomicNumber);
01713
01714 if(verboseLevel > 1)
01715 { G4cout << " After Tuning: " << G4endl;
01716 for(i=0; i<vecLen; i++) pv[i].Print(i);
01717 }
01718
01719
01720
01721 G4double tof = incidentTOF;
01722 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
01723 && (incidentKineticEnergy <= 0.2) )
01724 tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
01725
01726 for(i=0; i<vecLen; i++)
01727 {
01728 if(pv[i].getName() == "KaonZero" || pv[i].getName() == "AntiKaonZero")
01729 {
01730 pvmx[0] = pv[i];
01731 if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
01732 else pv[i].setDefinition("KaonZeroLong");
01733 pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
01734 }
01735 }
01736
01737 successful = true;
01738 delete [] pvmx;
01739 G4int testCurr=0;
01740 G4double totKin=0;
01741 for(testCurr=0; testCurr<vecLen; testCurr++)
01742 {
01743 totKin+=pv[testCurr].getKineticEnergy();
01744 if(totKin>incidentKineticEnergy*1.05)
01745 {
01746 vecLen = testCurr;
01747 break;
01748 }
01749 }
01750
01751 return;
01752 }
01753
01754 void
01755 G4HEInelastic::TuningOfHighEnergyCascading(G4HEVector pv[],
01756 G4int& vecLen,
01757 const G4HEVector& incidentParticle,
01758 const G4HEVector& targetParticle,
01759 G4double atomicWeight,
01760 G4double atomicNumber)
01761 {
01762 G4int i,j;
01763 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
01764 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
01765 G4double incidentCharge = incidentParticle.getCharge();
01766 G4double incidentMass = incidentParticle.getMass();
01767 G4double targetMass = targetParticle.getMass();
01768 G4int pionPlusCode = PionPlus.getCode();
01769 G4int pionMinusCode = PionMinus.getCode();
01770 G4int pionZeroCode = PionZero.getCode();
01771 G4int protonCode = Proton.getCode();
01772 G4int neutronCode = Neutron.getCode();
01773 G4HEVector *pvmx = new G4HEVector [10];
01774 G4double *reddec = new G4double [7];
01775
01776 if (incidentKineticEnergy > (25.+G4UniformRand()*75.) ) {
01777 G4double reden = -0.7 + 0.29*std::log10(incidentKineticEnergy);
01778
01779
01780 G4double redat = 0.722 - 0.278*std::log10(atomicWeight);
01781 G4double pmax = -200.;
01782 G4double pmapim = -200.;
01783 G4double pmapi0 = -200.;
01784 G4double pmapip = -200.;
01785 G4int ipmax = 0;
01786 G4int maxpim = 0;
01787 G4int maxpi0 = 0;
01788 G4int maxpip = 0;
01789 G4int iphmf;
01790 if ( (G4UniformRand() > (atomicWeight/100.-0.28))
01791 && (std::fabs(incidentCharge) > 0.) )
01792 {
01793 for (i=0; i < vecLen; i++)
01794 {
01795 iphmf = pv[i].getCode();
01796 G4double ppp = pv[i].Length();
01797 if ( ppp > pmax)
01798 {
01799 pmax = ppp; ipmax = i;
01800 }
01801 if (iphmf == pionPlusCode && ppp > pmapip )
01802 {
01803 pmapip = ppp; maxpip = i;
01804 }
01805 else if (iphmf == pionZeroCode && ppp > pmapi0)
01806 {
01807 pmapi0 = ppp; maxpi0 = i;
01808 }
01809 else if (iphmf == pionMinusCode && ppp > pmapim)
01810 {
01811 pmapim = ppp; maxpim = i;
01812 }
01813 }
01814 }
01815 if(verboseLevel > 1)
01816 {
01817 G4cout << "ipmax, pmax " << ipmax << " " << pmax << G4endl;
01818 G4cout << "maxpip,pmapip " << maxpip << " " << pmapip << G4endl;
01819 G4cout << "maxpi0,pmapi0 " << maxpi0 << " " << pmapi0 << G4endl;
01820 G4cout << "maxpim,pmapim " << maxpim << " " << pmapim << G4endl;
01821 }
01822
01823 if ( vecLen > 2)
01824 {
01825 for (i=2; i<vecLen; i++)
01826 {
01827 iphmf = pv[i].getCode();
01828 if ( ((iphmf==protonCode)||(iphmf==neutronCode)||(pv[i].getType()=="Nucleus"))
01829 && (pv[i].Length()<1.5)
01830 && ((G4UniformRand()<reden)||(G4UniformRand()<redat)))
01831 {
01832 pv[i].setMomentumAndUpdate( 0., 0., 0.);
01833 if(verboseLevel > 1)
01834 {
01835 G4cout << "zero Momentum for particle " << G4endl;
01836 pv[i].Print(i);
01837 }
01838 }
01839 }
01840 }
01841 if (maxpi0 == ipmax)
01842 {
01843 if (G4UniformRand() < pmapi0/incidentTotalMomentum)
01844 {
01845 if ((incidentCharge > 0.5) && (maxpip != 0))
01846 {
01847 G4ParticleMomentum mompi0 = pv[maxpi0].getMomentum();
01848 pv[maxpi0].setMomentumAndUpdate( pv[maxpip].getMomentum() );
01849 pv[maxpip].setMomentumAndUpdate( mompi0);
01850 if(verboseLevel > 1)
01851 {
01852 G4cout << " exchange Momentum for " << maxpi0 << " and " << maxpip << G4endl;
01853 }
01854 }
01855 else if ((incidentCharge < -0.5) && (maxpim != 0))
01856 {
01857 G4ParticleMomentum mompi0 = pv[maxpi0].getMomentum();
01858 pv[maxpi0].setMomentumAndUpdate( pv[maxpim].getMomentum() );
01859 pv[maxpim].setMomentumAndUpdate( mompi0);
01860 if(verboseLevel > 1)
01861 {
01862 G4cout << " exchange Momentum for " << maxpi0 << " and " << maxpip << G4endl;
01863 }
01864 }
01865 }
01866 }
01867 G4double bntot = - incidentParticle.getBaryonNumber() - atomicWeight;
01868 for (i=0; i<vecLen; i++) bntot += pv[i].getBaryonNumber();
01869 if(atomicWeight < 1.5) { bntot = 0.; }
01870 else { bntot = 1. + bntot/atomicWeight; }
01871 if(atomicWeight > (75.+G4UniformRand()*50.)) bntot = 0.;
01872 if(verboseLevel > 1)
01873 {
01874 G4cout << " Calculated Baryon- Number " << bntot << G4endl;
01875 }
01876
01877 j = 0;
01878 for (i=0; i<vecLen; i++)
01879 {
01880 G4double ppp = pv[i].Length();
01881 if ( ppp > 1.e-6)
01882 {
01883 iphmf = pv[i].getCode();
01884 if( (bntot > 0.3)
01885 && ((iphmf == protonCode) || (iphmf == neutronCode)
01886 || (pv[i].getType() == "Nucleus") )
01887 && (G4UniformRand() < 0.25)
01888 && (ppp < 1.2) )
01889 {
01890 if(verboseLevel > 1)
01891 {
01892 G4cout << " skip Baryon: " << G4endl;
01893 pv[i].Print(i);
01894 }
01895 continue;
01896
01897 }
01898 if (j != i)
01899 {
01900 pv[j] = pv[i];
01901 }
01902 j++;
01903 }
01904 }
01905 vecLen = j;
01906 }
01907
01908 pvmx[0] = incidentParticle;
01909 pvmx[1] = targetParticle;
01910 pvmx[8].setZero();
01911 pvmx[2].Add(pvmx[0], pvmx[1]);
01912 pvmx[3].Lor(pvmx[0], pvmx[2]);
01913 pvmx[4].Lor(pvmx[1], pvmx[2]);
01914
01915 if (verboseLevel > 1) {
01916 pvmx[0].Print(0);
01917 incidentParticle.Print(0);
01918 pvmx[1].Print(1);
01919 targetParticle.Print(1);
01920 pvmx[2].Print(2);
01921 pvmx[3].Print(3);
01922 pvmx[4].Print(4);
01923 }
01924
01925
01926
01927
01928
01929
01930
01931 G4int ledpar = -1;
01932 G4double redpar = 0.;
01933 G4int incidentS = incidentParticle.getStrangenessNumber();
01934 if (incidentParticle.getName() == "KaonZeroShort" ||
01935 incidentParticle.getName() == "KaonZeroLong") {
01936 if(G4UniformRand() < 0.5) {
01937 incidentS = 1;
01938 } else {
01939 incidentS = -1;
01940 }
01941 }
01942
01943 G4int incidentB = incidentParticle.getBaryonNumber();
01944
01945 for (i=0; i<vecLen; i++) {
01946 G4int iphmf = pv[i].getCode();
01947 G4double ppp = pv[i].Length();
01948
01949 if (ppp > 1.e-3) {
01950 pvmx[5].Lor( pv[i], pvmx[2] );
01951 G4double cost = pvmx[3].CosAng( pvmx[5] );
01952
01953
01954
01955
01956
01957 G4int particleS = pv[i].getStrangenessNumber();
01958
01959 if (pv[i].getName() == "KaonZeroShort" ||
01960 pv[i].getName() == "KaonZeroLong") {
01961 if (G4UniformRand() < 0.5) {
01962 particleS = 1;
01963 } else {
01964 particleS = -1;
01965 }
01966 }
01967 G4int particleB = pv[i].getBaryonNumber();
01968 G4double hfmass;
01969 if (cost > 0.) {
01970 reddec[0] = std::fabs( incidentMass - pv[i].getMass() );
01971 reddec[1] = std::fabs( incidentCharge - pv[i].getCharge());
01972 reddec[2] = std::fabs( G4double(incidentS - particleS) );
01973 reddec[3] = std::fabs( G4double(incidentB - particleB) );
01974 hfmass = incidentMass;
01975
01976 } else {
01977 reddec[0] = std::fabs( targetMass - pv[i].getMass() );
01978 reddec[1] = std::fabs( atomicNumber/atomicWeight - pv[i].getCharge());
01979 reddec[2] = std::fabs( G4double(particleS) );
01980 reddec[3] = std::fabs( 1. - particleB );
01981 hfmass = targetMass;
01982 }
01983
01984 reddec[5] = reddec[0]+reddec[1]+reddec[2]+reddec[3];
01985 G4double sbqwgt = reddec[5];
01986 if (hfmass < 0.2) {
01987 sbqwgt = (sbqwgt-2.5)*0.10;
01988 if(pv[i].getCode() == pionZeroCode) sbqwgt = 0.15;
01989 } else if (hfmass < 0.6) {
01990 sbqwgt = (sbqwgt-3.0)*0.10;
01991 } else {
01992 sbqwgt = (sbqwgt-2.0)*0.10;
01993 if(pv[i].getCode() == pionZeroCode) sbqwgt = 0.15;
01994 }
01995
01996 ppp = pvmx[5].Length();
01997
01998
01999
02000
02001 if (sbqwgt > 0. && ppp > 1.e-6) {
02002 G4double pthmf = ppp*std::sqrt(1.-cost*cost);
02003 G4double plhmf = ppp*cost*(1.-sbqwgt);
02004 pvmx[7].Cross( pvmx[3], pvmx[5] );
02005 pvmx[7].Cross( pvmx[7], pvmx[3] );
02006
02007 if (pvmx[3].Length() > 0.)
02008 pvmx[6].SmulAndUpdate( pvmx[3], plhmf/pvmx[3].Length() );
02009 else if(verboseLevel > 1)
02010 G4cout << "NaNQ in Tuning of High Energy Hadronic Interactions" << G4endl;
02011
02012 if (pvmx[7].Length() > 0.)
02013 pvmx[7].SmulAndUpdate( pvmx[7], pthmf/pvmx[7].Length() );
02014 else if(verboseLevel > 1)
02015 G4cout << "NaNQ in Tuning of High Energy Hadronic Interactions" << G4endl;
02016
02017 pvmx[5].Add3(pvmx[6], pvmx[7] );
02018 pvmx[5].setEnergy( std::sqrt(sqr(pvmx[5].Length()) + sqr(pv[i].getMass())));
02019 pv[i].Lor( pvmx[5], pvmx[4] );
02020 if (verboseLevel > 1) {
02021 G4cout << " Particle Momentum changed to: " << G4endl;
02022 pv[i].Print(i);
02023 }
02024 }
02025
02026
02027
02028
02029
02030 G4int ss = -3;
02031 if (incidentS != 0) ss = 0;
02032 if (iphmf != pionZeroCode && pv[i].getSide() > ss) {
02033 pvmx[7].Sub3( incidentParticle, pv[i] );
02034 reddec[4] = pvmx[7].Length()/incidentTotalMomentum;
02035 reddec[6] = reddec[4]*2./3. + reddec[5]/12.;
02036 reddec[6] = Amax(0., 1. - reddec[6]);
02037 if ( (reddec[5] <= 3.75) && (reddec[6] > redpar) ) {
02038 ledpar = i;
02039 redpar = reddec[6];
02040 }
02041 }
02042 }
02043 pvmx[8].Add3(pvmx[8], pv[i] );
02044 }
02045
02046 if(false) if (ledpar >= 0)
02047 {
02048 if(verboseLevel > 1)
02049 {
02050 G4cout << " Leading Particle found : " << ledpar << G4endl;
02051 pv[ledpar].Print(ledpar);
02052 pvmx[8].Print(-2);
02053 incidentParticle.Print(-1);
02054 }
02055 pvmx[4].Sub3(incidentParticle,pvmx[8]);
02056 pvmx[5].Smul(incidentParticle, incidentParticle.CosAng(pvmx[4])
02057 *pvmx[4].Length()/incidentParticle.Length());
02058 pv[ledpar].Add3(pv[ledpar],pvmx[5]);
02059
02060 pv[ledpar].SmulAndUpdate( pv[ledpar], 1.);
02061 if(verboseLevel > 1)
02062 {
02063 pv[ledpar].Print(ledpar);
02064 }
02065 }
02066
02067 if (conserveEnergy) {
02068 G4double ekinhf = 0.;
02069 for (i=0; i<vecLen; i++) {
02070 ekinhf += pv[i].getKineticEnergy();
02071 if(pv[i].getMass() < 0.7) ekinhf += pv[i].getMass();
02072 }
02073 if(incidentParticle.getMass() < 0.7) ekinhf -= incidentParticle.getMass();
02074
02075 if(ledpar < 0) {
02076 ekinhf = incidentParticle.getKineticEnergy()/ekinhf;
02077 for (i=0; i<vecLen; i++)
02078 pv[i].setKineticEnergyAndUpdate(ekinhf*pv[i].getKineticEnergy());
02079
02080 } else {
02081
02082
02083 ekinhf = incidentParticle.getKineticEnergy() - ekinhf;
02084 ekinhf += pv[ledpar].getKineticEnergy();
02085 if(ekinhf < 0.) ekinhf = 0.;
02086 pv[ledpar].setKineticEnergyAndUpdate(ekinhf);
02087 }
02088 }
02089
02090 delete [] reddec;
02091 delete [] pvmx;
02092
02093 return;
02094 }
02095
02096 void
02097 G4HEInelastic::HighEnergyClusterProduction(G4bool& successful,
02098 G4HEVector pv[],
02099 G4int& vecLen,
02100 G4double& excitationEnergyGNP,
02101 G4double& excitationEnergyDTA,
02102 const G4HEVector& incidentParticle,
02103 const G4HEVector& targetParticle,
02104 G4double atomicWeight,
02105 G4double atomicNumber)
02106 {
02107
02108
02109
02110
02111
02112
02113
02114 G4int protonCode = Proton.getCode();
02115 G4double protonMass = Proton.getMass();
02116 G4int neutronCode = Neutron.getCode();
02117 G4double kaonPlusMass = KaonPlus.getMass();
02118 G4int pionPlusCode = PionPlus.getCode();
02119 G4int pionZeroCode = PionZero.getCode();
02120 G4int pionMinusCode = PionMinus.getCode();
02121 G4String mesonType = PionPlus.getType();
02122 G4String baryonType = Proton.getType();
02123 G4String antiBaryonType = AntiProton.getType();
02124
02125 G4double targetMass = targetParticle.getMass();
02126
02127 G4int incidentCode = incidentParticle.getCode();
02128 G4double incidentMass = incidentParticle.getMass();
02129 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
02130 G4double incidentEnergy = incidentParticle.getEnergy();
02131 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
02132 G4String incidentType = incidentParticle.getType();
02133
02134 G4double incidentTOF = 0.;
02135
02136
02137 G4int i, j;
02138
02139 if (verboseLevel > 1)
02140 G4cout << " G4HEInelastic::HighEnergyClusterProduction " << G4endl;
02141
02142 successful = false;
02143 if (incidentTotalMomentum < 25. + G4UniformRand()*25.) return;
02144
02145 G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
02146 +2.*targetMass*incidentEnergy);
02147
02148 G4HEVector pvI = incidentParticle;
02149 pvI.setSide(1);
02150
02151 G4HEVector pvT = targetParticle;
02152 pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
02153 pvT.setSide( -1 );
02154 pvT.setTOF( -1.);
02155
02156
02157
02158
02159 G4int targ = 0;
02160 G4int ifor = 0;
02161 G4int iback = 0;
02162 G4int pvCode;
02163 G4double pvMass, pvEnergy;
02164
02165 pv[0].setSide(1);
02166 pv[1].setSide(-1);
02167 for (i = 0; i < vecLen; i++) {
02168 if (i > 1) {
02169 if (G4UniformRand() < 0.5) {
02170 pv[i].setSide( 1 );
02171 if (++ifor > 18) {
02172 pv[i].setSide(-1);
02173 ifor--;
02174 iback++;
02175 }
02176 } else {
02177 pv[i].setSide( -1 );
02178 if (++iback > 18) {
02179 pv[i].setSide(1);
02180 ifor++;
02181 iback--;
02182 }
02183 }
02184 }
02185
02186 pvCode = pv[i].getCode();
02187
02188 if ( ( (incidentCode == protonCode) || (incidentCode == neutronCode)
02189 || (incidentType == mesonType) )
02190 && ( (pvCode == pionPlusCode) || (pvCode == pionMinusCode) )
02191 && ( (G4UniformRand() < (10.-incidentTotalMomentum)/6.) )
02192 && ( (G4UniformRand() < atomicWeight/300.) ) ) {
02193 if (G4UniformRand() > atomicNumber/atomicWeight)
02194 pv[i].setDefinition("Neutron");
02195 else
02196 pv[i].setDefinition("Proton");
02197 targ++;
02198 }
02199 pv[i].setTOF( incidentTOF );
02200 }
02201
02202 G4double tb = 2. * iback;
02203 if (centerOfMassEnergy < (2+G4UniformRand())) tb = (2.*iback + vecLen)/2.;
02204
02205 G4double nucsup[] = {1.0, 0.7, 0.5, 0.4, 0.35, 0.3};
02206 G4double psup[] = {3. , 6. , 20., 50., 100.,1000.};
02207 G4double ss = centerOfMassEnergy*centerOfMassEnergy;
02208
02209 G4double xtarg = Amax(0.01, Amin(0.50, 0.312+0.2*std::log(std::log(ss))+std::pow(ss,1.5)/6000.)
02210 * (std::pow(atomicWeight,0.33)-1.) * tb);
02211 G4int momentumBin = 0;
02212 while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin])) momentumBin++;
02213 momentumBin = Imin(5, momentumBin);
02214 G4double xpnhmf = Amax(0.01,xtarg*nucsup[momentumBin]);
02215 G4double xshhmf = Amax(0.01,xtarg-xpnhmf);
02216 G4double rshhmf = 0.25*xshhmf;
02217 G4double rpnhmf = 0.81*xpnhmf;
02218 G4double xhmf;
02219 G4int nshhmf, npnhmf;
02220 if (rshhmf > 1.1)
02221 {
02222 rshhmf = xshhmf/(rshhmf-1.);
02223 if (rshhmf <= 20.)
02224 xhmf = GammaRand( rshhmf );
02225 else
02226 xhmf = Erlang( G4int(rshhmf+0.5) );
02227 xshhmf *= xhmf/rshhmf;
02228 }
02229 nshhmf = Poisson( xshhmf );
02230 if (rpnhmf > 1.1)
02231 {
02232 rpnhmf = xpnhmf/(rpnhmf -1.);
02233 if (rpnhmf <= 20.)
02234 xhmf = GammaRand( rpnhmf );
02235 else
02236 xhmf = Erlang( G4int(rpnhmf+0.5) );
02237 xpnhmf *= xhmf/rpnhmf;
02238 }
02239 npnhmf = Poisson( xpnhmf );
02240
02241 while (npnhmf > 0)
02242 {
02243 if ( G4UniformRand() > (1. - atomicNumber/atomicWeight))
02244 pv[vecLen].setDefinition( "Proton" );
02245 else
02246 pv[vecLen].setDefinition( "Neutron" );
02247 targ++;
02248 pv[vecLen].setSide( -2 );
02249 pv[vecLen].setFlag( true );
02250 pv[vecLen].setTOF( incidentTOF );
02251 vecLen++;
02252 npnhmf--;
02253 }
02254 while (nshhmf > 0)
02255 {
02256 G4double ran = G4UniformRand();
02257 if (ran < 0.333333 )
02258 pv[vecLen].setDefinition( "PionPlus" );
02259 else if (ran < 0.6667)
02260 pv[vecLen].setDefinition( "PionZero" );
02261 else
02262 pv[vecLen].setDefinition( "PionMinus" );
02263 pv[vecLen].setSide( -2 );
02264 pv[vecLen].setFlag( true );
02265 pv[vecLen].setTOF( incidentTOF );
02266 vecLen++;
02267 nshhmf--;
02268 }
02269
02270
02271
02272
02273
02274
02275 G4int lead = 0;
02276 G4HEVector leadParticle;
02277 if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
02278 && (incidentCode != neutronCode) )
02279 {
02280 G4double pMass = pv[0].getMass();
02281 G4int pCode = pv[0].getCode();
02282 if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
02283 && (pCode != neutronCode) )
02284 {
02285 lead = pCode;
02286 leadParticle = pv[0];
02287 }
02288 else
02289 {
02290 pMass = pv[1].getMass();
02291 pCode = pv[1].getCode();
02292 if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
02293 && (pCode != neutronCode) )
02294 {
02295 lead = pCode;
02296 leadParticle = pv[1];
02297 }
02298 }
02299 }
02300
02301 if (verboseLevel > 1)
02302 { G4cout << " pv Vector after initialization " << vecLen << G4endl;
02303 pvI.Print(-1);
02304 pvT.Print(-1);
02305 for (i=0; i < vecLen ; i++) pv[i].Print(i);
02306 }
02307
02308 G4double tavai = 0.;
02309 for(i=0;i<vecLen;i++) if(pv[i].getSide() != -2) tavai += pv[i].getMass();
02310
02311 while (tavai > centerOfMassEnergy)
02312 {
02313 for (i=vecLen-1; i >= 0; i--)
02314 {
02315 if (pv[i].getSide() != -2)
02316 {
02317 tavai -= pv[i].getMass();
02318 if( i != vecLen-1)
02319 {
02320 for (j=i; j < vecLen; j++)
02321 {
02322 pv[j] = pv[j+1];
02323 }
02324 }
02325 if ( --vecLen < 2)
02326 {
02327 successful = false;
02328 return;
02329 }
02330 break;
02331 }
02332 }
02333 }
02334
02335
02336
02337
02338
02339
02340 G4double rmc0 = 0., rmd0 = 0., rme0 = 0.;
02341 G4int ntc = 0, ntd = 0, nte = 0;
02342
02343 for (i=0; i < vecLen; i++)
02344 {
02345 if(pv[i].getSide() > 0)
02346 {
02347 if(ntc < 17)
02348 {
02349 rmc0 += pv[i].getMass();
02350 ntc++;
02351 }
02352 else
02353 {
02354 if(ntd < 17)
02355 {
02356 pv[i].setSide(-1);
02357 rmd0 += pv[i].getMass();
02358 ntd++;
02359 }
02360 else
02361 {
02362 pv[i].setSide(-2);
02363 rme0 += pv[i].getMass();
02364 nte++;
02365 }
02366 }
02367 }
02368 else if (pv[i].getSide() == -1)
02369 {
02370 if(ntd < 17)
02371 {
02372 rmd0 += pv[i].getMass();
02373 ntd++;
02374 }
02375 else
02376 {
02377 pv[i].setSide(-2);
02378 rme0 += pv[i].getMass();
02379 nte++;
02380 }
02381 }
02382 else
02383 {
02384 rme0 += pv[i].getMass();
02385 nte++;
02386 }
02387 }
02388
02389 G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
02390 G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
02391
02392 G4double rmc = rmc0, rmd = rmd0, rme = rme0;
02393 G4int ntc1 = Imin(4,ntc-1);
02394 G4int ntd1 = Imin(4,ntd-1);
02395 G4int nte1 = Imin(4,nte-1);
02396 if (ntc > 1) rmc = rmc0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntc1])/gpar[ntc1];
02397 if (ntd > 1) rmd = rmd0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntd1])/gpar[ntd1];
02398 if (nte > 1) rme = rme0 + std::pow(-std::log(1.-G4UniformRand()),cpar[nte1])/gpar[nte1];
02399 while( (rmc+rmd) > centerOfMassEnergy)
02400 {
02401 if ((rmc == rmc0) && (rmd == rmd0))
02402 {
02403 rmd *= 0.999*centerOfMassEnergy/(rmc+rmd);
02404 rmc *= 0.999*centerOfMassEnergy/(rmc+rmd);
02405 }
02406 else
02407 {
02408 rmc = 0.1*rmc0 + 0.9*rmc;
02409 rmd = 0.1*rmd0 + 0.9*rmd;
02410 }
02411 }
02412 if(verboseLevel > 1)
02413 G4cout << " Cluster Masses: " << ntc << " " << rmc << " " << ntd
02414 << " " << rmd << " " << nte << " " << rme << G4endl;
02415
02416 G4HEVector* pvmx = new G4HEVector[11];
02417
02418 pvmx[1].setMass( incidentMass);
02419 pvmx[1].setMomentumAndUpdate( 0., 0., incidentTotalMomentum);
02420 pvmx[2].setMass( targetMass);
02421 pvmx[2].setMomentumAndUpdate( 0., 0., 0.);
02422 pvmx[0].Add( pvmx[1], pvmx[2] );
02423 pvmx[1].Lor( pvmx[1], pvmx[0] );
02424 pvmx[2].Lor( pvmx[2], pvmx[0] );
02425
02426 G4double pf = std::sqrt(Amax(0.0001, sqr(sqr(centerOfMassEnergy) + rmd*rmd -rmc*rmc)
02427 - 4*sqr(centerOfMassEnergy)*rmd*rmd))/(2.*centerOfMassEnergy);
02428 pvmx[3].setMass( rmc );
02429 pvmx[4].setMass( rmd );
02430 pvmx[3].setEnergy( std::sqrt(pf*pf + rmc*rmc) );
02431 pvmx[4].setEnergy( std::sqrt(pf*pf + rmd*rmd) );
02432
02433 G4double tvalue = -DBL_MAX;
02434 G4double bvalue = Amax(0.01, 4.0 + 1.6*std::log(incidentTotalMomentum));
02435 if (bvalue != 0.0) tvalue = std::log(G4UniformRand())/bvalue;
02436 G4double pin = pvmx[1].Length();
02437 G4double tacmin = sqr( pvmx[1].getEnergy() - pvmx[3].getEnergy()) - sqr( pin - pf);
02438 G4double ctet = Amax(-1., Amin(1., 1.+2.*(tvalue-tacmin)/Amax(1.e-10, 4.*pin*pf)));
02439 G4double stet = std::sqrt(Amax(0., 1.0 - ctet*ctet));
02440 G4double phi = twopi * G4UniformRand();
02441 pvmx[3].setMomentum( pf * stet * std::sin(phi),
02442 pf * stet * std::cos(phi),
02443 pf * ctet );
02444 pvmx[4].Smul( pvmx[3], -1.);
02445
02446 if (nte > 0)
02447 {
02448 G4double ekit1 = 0.04;
02449 G4double ekit2 = 0.6;
02450 G4double gaval = 1.2;
02451 if (incidentKineticEnergy <= 5.)
02452 {
02453 ekit1 *= sqr(incidentKineticEnergy)/25.;
02454 ekit2 *= sqr(incidentKineticEnergy)/25.;
02455 }
02456 G4double avalue = (1.-gaval)/(std::pow(ekit2, 1.-gaval)-std::pow(ekit1, 1.-gaval));
02457 for (i=0; i < vecLen; i++)
02458 {
02459 if (pv[i].getSide() == -2)
02460 {
02461 G4double ekit = std::pow(G4UniformRand()*(1.-gaval)/avalue +std::pow(ekit1, 1.-gaval),
02462 1./(1.-gaval));
02463 pv[i].setKineticEnergyAndUpdate( ekit );
02464 ctet = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
02465 stet = std::sqrt( Amax( 0.0, 1. - ctet*ctet ));
02466 phi = G4UniformRand()*twopi;
02467 G4double pp = pv[i].Length();
02468 pv[i].setMomentum( pp * stet * std::sin(phi),
02469 pp * stet * std::cos(phi),
02470 pp * ctet );
02471 pv[i].Lor( pv[i], pvmx[0] );
02472 }
02473 }
02474 }
02475
02476
02477 pvmx[5].SmulAndUpdate( pvmx[3], -1.);
02478 pvmx[6].SmulAndUpdate( pvmx[4], -1.);
02479
02480 if (verboseLevel > 1) {
02481 G4cout << " General vectors before Phase space Generation " << G4endl;
02482 for (i=0; i<7; i++) pvmx[i].Print(i);
02483 }
02484
02485 G4HEVector* tempV = new G4HEVector[18];
02486 G4bool constantCrossSection = true;
02487 G4double wgt;
02488 G4int npg;
02489
02490 if (ntc > 1)
02491 {
02492 npg = 0;
02493 for (i=0; i < vecLen; i++)
02494 {
02495 if (pv[i].getSide() > 0)
02496 {
02497 tempV[npg++] = pv[i];
02498 }
02499 }
02500 wgt = NBodyPhaseSpace( pvmx[3].getMass(), constantCrossSection, tempV, npg);
02501
02502 npg = 0;
02503 for (i=0; i < vecLen; i++)
02504 {
02505 if (pv[i].getSide() > 0)
02506 {
02507 pv[i].setMomentum( tempV[npg++].getMomentum());
02508 pv[i].SmulAndUpdate( pv[i], 1. );
02509 pv[i].Lor( pv[i], pvmx[5] );
02510 }
02511 }
02512 }
02513 else if(ntc == 1)
02514 {
02515 for(i=0; i<vecLen; i++)
02516 {
02517 if(pv[i].getSide() > 0) pv[i].setMomentumAndUpdate(pvmx[3].getMomentum());
02518 }
02519 }
02520 else
02521 {
02522 }
02523
02524 if (ntd > 1)
02525 {
02526 npg = 0;
02527 for (i=0; i < vecLen; i++)
02528 {
02529 if (pv[i].getSide() == -1)
02530 {
02531 tempV[npg++] = pv[i];
02532 }
02533 }
02534 wgt = NBodyPhaseSpace( pvmx[4].getMass(), constantCrossSection, tempV, npg);
02535
02536 npg = 0;
02537 for (i=0; i < vecLen; i++)
02538 {
02539 if (pv[i].getSide() == -1)
02540 {
02541 pv[i].setMomentum( tempV[npg++].getMomentum());
02542 pv[i].SmulAndUpdate( pv[i], 1.);
02543 pv[i].Lor( pv[i], pvmx[6] );
02544 }
02545 }
02546 }
02547 else if(ntd == 1)
02548 {
02549 for(i=0; i<vecLen; i++)
02550 {
02551 if(pv[i].getSide() == -1) pv[i].setMomentumAndUpdate(pvmx[4].getMomentum());
02552 }
02553 }
02554 else
02555 {
02556 }
02557
02558 if(verboseLevel > 1)
02559 {
02560 G4cout << " Vectors after PhaseSpace generation " << G4endl;
02561 for(i=0; i<vecLen; i++) pv[i].Print(i);
02562 }
02563
02564
02565
02566 targ = 0;
02567 for( i=0; i < vecLen; i++ )
02568 {
02569 if( pv[i].getType() == baryonType )targ++;
02570 if( pv[i].getType() == antiBaryonType )targ--;
02571 pv[i].Lor( pv[i], pvmx[2] );
02572 }
02573 if (targ<1) targ = 1;
02574
02575 if(verboseLevel > 1) {
02576 G4cout << " Transformation in Lab- System " << G4endl;
02577 for(i=0; i<vecLen; i++) pv[i].Print(i);
02578 }
02579
02580 G4bool dum(0);
02581 G4double ekin, teta;
02582
02583 if( lead )
02584 {
02585 for( i=0; i<vecLen; i++ )
02586 {
02587 if( pv[i].getCode() == lead )
02588 {
02589 dum = false;
02590 break;
02591 }
02592 }
02593 if( dum )
02594 {
02595 i = 0;
02596
02597 if( ( (leadParticle.getType() == baryonType ||
02598 leadParticle.getType() == antiBaryonType)
02599 && (pv[1].getType() == baryonType ||
02600 pv[1].getType() == antiBaryonType))
02601 || ( (leadParticle.getType() == mesonType)
02602 && (pv[1].getType() == mesonType)))
02603 {
02604 i = 1;
02605 }
02606
02607 ekin = pv[i].getKineticEnergy();
02608 pv[i] = leadParticle;
02609 if( pv[i].getFlag() )
02610 pv[i].setTOF( -1.0 );
02611 else
02612 pv[i].setTOF( 1.0 );
02613 pv[i].setKineticEnergyAndUpdate( ekin );
02614 }
02615 }
02616
02617 pvmx[4].setMass( incidentMass);
02618 pvmx[4].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
02619
02620 G4double ekin0 = pvmx[4].getKineticEnergy();
02621
02622 pvmx[5].setMass ( protonMass * targ);
02623 pvmx[5].setEnergy( protonMass * targ);
02624 pvmx[5].setKineticEnergy(0.);
02625 pvmx[5].setMomentum( 0.0, 0.0, 0.0 );
02626
02627 ekin = pvmx[4].getEnergy() + pvmx[5].getEnergy();
02628
02629 pvmx[6].Add( pvmx[4], pvmx[5] );
02630 pvmx[4].Lor( pvmx[4], pvmx[6] );
02631 pvmx[5].Lor( pvmx[5], pvmx[6] );
02632
02633 G4double tecm = pvmx[4].getEnergy() + pvmx[5].getEnergy();
02634
02635 pvmx[8].setZero();
02636
02637 G4double ekin1 = 0.0;
02638
02639 for( i=0; i < vecLen; i++ )
02640 {
02641 pvmx[8].Add( pvmx[8], pv[i] );
02642 ekin1 += pv[i].getKineticEnergy();
02643 ekin -= pv[i].getMass();
02644 }
02645
02646 if( vecLen > 1 && vecLen < 19 )
02647 {
02648 constantCrossSection = true;
02649 G4HEVector pw[18];
02650 for(i=0;i<vecLen;i++) pw[i] = pv[i];
02651 wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
02652 ekin = 0.0;
02653 for( i=0; i < vecLen; i++ )
02654 {
02655 pvmx[7].setMass( pw[i].getMass());
02656 pvmx[7].setMomentum( pw[i].getMomentum() );
02657 pvmx[7].SmulAndUpdate( pvmx[7], 1.);
02658 pvmx[7].Lor( pvmx[7], pvmx[5] );
02659 ekin += pvmx[7].getKineticEnergy();
02660 }
02661 teta = pvmx[8].Ang( pvmx[4] );
02662 if (verboseLevel > 1)
02663 G4cout << " vecLen > 1 && vecLen < 19 " << teta << " "
02664 << ekin0 << " " << ekin1 << " " << ekin << G4endl;
02665 }
02666
02667 if( ekin1 != 0.0 )
02668 {
02669 pvmx[7].setZero();
02670 wgt = ekin/ekin1;
02671 ekin1 = 0.;
02672 for( i=0; i < vecLen; i++ )
02673 {
02674 pvMass = pv[i].getMass();
02675 ekin = pv[i].getKineticEnergy() * wgt;
02676 pv[i].setKineticEnergyAndUpdate( ekin );
02677 ekin1 += ekin;
02678 pvmx[7].Add( pvmx[7], pv[i] );
02679 }
02680 teta = pvmx[7].Ang( pvmx[4] );
02681 if (verboseLevel > 1)
02682 G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
02683 << ekin1 << G4endl;
02684 }
02685
02686 if(verboseLevel > 1)
02687 {
02688 G4cout << " After energy- correction " << G4endl;
02689 for(i=0; i<vecLen; i++) pv[i].Print(i);
02690 }
02691
02692
02693
02694 G4double ry = G4UniformRand();
02695 G4double rz = G4UniformRand();
02696 G4double rx = twopi*rz;
02697 G4double a1 = std::sqrt(-2.0*std::log(ry));
02698 G4double rantarg1 = a1*std::cos(rx)*0.02*targ/G4double(vecLen);
02699 G4double rantarg2 = a1*std::sin(rx)*0.02*targ/G4double(vecLen);
02700
02701 for (i = 0; i < vecLen; i++)
02702 pv[i].setMomentum(pv[i].getMomentum().x()+rantarg1,
02703 pv[i].getMomentum().y()+rantarg2);
02704
02705 if (verboseLevel > 1) {
02706 pvmx[7].setZero();
02707 for (i = 0; i < vecLen; i++) pvmx[7].Add( pvmx[7], pv[i] );
02708 teta = pvmx[7].Ang( pvmx[4] );
02709 G4cout << " After smearing " << teta << G4endl;
02710 }
02711
02712
02713
02714
02715
02716
02717
02718
02719 G4double dekin = 0.0;
02720 G4int npions = 0;
02721 G4double ek1 = 0.0;
02722 G4double alekw, xxh;
02723 G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
02724 G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.0};
02725 G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
02726
02727 for (i = 0; i < vecLen; i++) {
02728 pv[i].Defs1( pv[i], pvI );
02729 if (atomicWeight > 1.5) {
02730 ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa*( 1. + 0.5*normal()));
02731 alekw = std::log( incidentKineticEnergy );
02732 xxh = 1.;
02733 if (incidentCode == pionPlusCode || incidentCode == pionMinusCode) {
02734 if (pv[i].getCode() == pionZeroCode) {
02735 if (G4UniformRand() < std::log(atomicWeight)) {
02736 if (alekw > alem[0]) {
02737 G4int jmax = 1;
02738 for (j = 1; j < 8; j++) {
02739 if (alekw < alem[j]) {
02740 jmax = j;
02741 break;
02742 }
02743 }
02744 xxh = (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alekw
02745 + val0[jmax-1] - (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alem[jmax-1];
02746 xxh = 1. - xxh;
02747 }
02748 }
02749 }
02750 }
02751 dekin += ekin*(1.-xxh);
02752 ekin *= xxh;
02753 pv[i].setKineticEnergyAndUpdate( ekin );
02754 pvCode = pv[i].getCode();
02755 if ((pvCode == pionPlusCode) ||
02756 (pvCode == pionMinusCode) ||
02757 (pvCode == pionZeroCode)) {
02758 npions += 1;
02759 ek1 += ekin;
02760 }
02761 }
02762 }
02763
02764 if( (ek1 > 0.0) && (npions > 0) )
02765 {
02766 dekin = 1.+dekin/ek1;
02767 for (i = 0; i < vecLen; i++)
02768 {
02769 pvCode = pv[i].getCode();
02770 if((pvCode == pionPlusCode) || (pvCode == pionMinusCode) || (pvCode == pionZeroCode))
02771 {
02772 ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
02773 pv[i].setKineticEnergyAndUpdate( ekin );
02774 }
02775 }
02776 }
02777 if (verboseLevel > 1)
02778 { G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
02779 for (i=0; i<vecLen; i++) pv[i].Print(i);
02780 }
02781
02782
02783
02784
02785
02786 if (verboseLevel > 1)
02787 G4cout << " Evaporation " << atomicWeight << " " << excitationEnergyGNP
02788 << " " << excitationEnergyDTA << G4endl;
02789
02790 G4double sprob = 0.;
02791 if (incidentKineticEnergy > 5. )
02792
02793 sprob = Amin(1., 0.000314*atomicWeight*std::log(incidentKineticEnergy-4.));
02794 if( atomicWeight > 1.5 && G4UniformRand() > sprob)
02795 {
02796
02797 G4double cost, sint, ekin2, ran, pp, eka;
02798 G4int spall(0), nbl(0);
02799
02800
02801
02802 if( excitationEnergyGNP >= 0.001 )
02803 {
02804
02805
02806
02807 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyGNP/
02808 (excitationEnergyGNP+excitationEnergyDTA));
02809 if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
02810 if (verboseLevel > 1)
02811 G4cout << " evaporation " << targ << " " << nbl << " "
02812 << sprob << G4endl;
02813 spall = targ;
02814 if( nbl > 0)
02815 {
02816 ekin = excitationEnergyGNP/nbl;
02817 ekin2 = 0.0;
02818 for( i=0; i<nbl; i++ )
02819 {
02820 if( G4UniformRand() < sprob ) continue;
02821 if( ekin2 > excitationEnergyGNP) break;
02822 ran = G4UniformRand();
02823 ekin1 = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
02824 if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
02825 ekin2 += ekin1;
02826 if( ekin2 > excitationEnergyGNP)
02827 ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
02828 if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
02829 pv[vecLen].setDefinition( "Proton");
02830 else
02831 pv[vecLen].setDefinition( "Neutron" );
02832 spall++;
02833 cost = G4UniformRand() * 2.0 - 1.0;
02834 sint = std::sqrt(std::fabs(1.0-cost*cost));
02835 phi = twopi * G4UniformRand();
02836 pv[vecLen].setFlag( true );
02837 pv[vecLen].setSide( -4 );
02838 pvMass = pv[vecLen].getMass();
02839 pv[vecLen].setTOF( 1.0 );
02840 pvEnergy = ekin1 + pvMass;
02841 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
02842 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
02843 pp*sint*std::cos(phi),
02844 pp*cost );
02845 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
02846 vecLen++;
02847 }
02848 if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
02849 {
02850 G4int ika, kk = 0;
02851 eka = incidentKineticEnergy;
02852 if( eka > 1.0 )eka *= eka;
02853 eka = Amax( 0.1, eka );
02854 ika = G4int(3.6*std::exp((atomicNumber*atomicNumber
02855 /atomicWeight-35.56)/6.45)/eka);
02856 if( ika > 0 )
02857 {
02858 for( i=(vecLen-1); i>=0; i-- )
02859 {
02860 if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
02861 {
02862 G4HEVector pTemp = pv[i];
02863 pv[i].setDefinition( "Neutron" );
02864 pv[i].setMomentumAndUpdate(pTemp.getMomentum());
02865 if (verboseLevel > 1) pv[i].Print(i);
02866 if( ++kk > ika ) break;
02867 }
02868 }
02869 }
02870 }
02871 }
02872 }
02873
02874
02875
02876
02877 if( excitationEnergyDTA >= 0.001 )
02878 {
02879 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyDTA
02880 /(excitationEnergyGNP+excitationEnergyDTA));
02881
02882
02883
02884 if( nbl > 0 )
02885 {
02886 ekin = excitationEnergyDTA/nbl;
02887 ekin2 = 0.0;
02888 for( i=0; i<nbl; i++ )
02889 {
02890 if( G4UniformRand() < sprob ) continue;
02891 if( ekin2 > excitationEnergyDTA) break;
02892 ran = G4UniformRand();
02893 ekin1 = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
02894 if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
02895 ekin2 += ekin1;
02896 if( ekin2 > excitationEnergyDTA )
02897 ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
02898 cost = G4UniformRand()*2.0 - 1.0;
02899 sint = std::sqrt(std::fabs(1.0-cost*cost));
02900 phi = twopi*G4UniformRand();
02901 ran = G4UniformRand();
02902 if( ran <= 0.60 )
02903 pv[vecLen].setDefinition( "Deuteron");
02904 else if (ran <= 0.90)
02905 pv[vecLen].setDefinition( "Triton" );
02906 else
02907 pv[vecLen].setDefinition( "Alpha" );
02908 spall += (int)(pv[vecLen].getMass() * 1.066);
02909 if( spall > atomicWeight ) break;
02910 pv[vecLen].setFlag( true );
02911 pv[vecLen].setSide( -4 );
02912 pvMass = pv[vecLen].getMass();
02913 pv[vecLen].setTOF( 1.0 );
02914 pvEnergy = pvMass + ekin1;
02915 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
02916 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
02917 pp*sint*std::cos(phi),
02918 pp*cost );
02919 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
02920 vecLen++;
02921 }
02922 }
02923 }
02924 }
02925 if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
02926 {
02927 for( i=0; i<vecLen; i++ )
02928 {
02929 G4double etb = pv[i].getKineticEnergy();
02930 if( etb >= incidentKineticEnergy )
02931 pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
02932 }
02933 }
02934
02935 TuningOfHighEnergyCascading( pv, vecLen,
02936 incidentParticle, targetParticle,
02937 atomicWeight, atomicNumber);
02938
02939
02940
02941 G4double tof = incidentTOF;
02942 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
02943 && (incidentKineticEnergy <= 0.2) )
02944 tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
02945 for ( i=0; i < vecLen; i++)
02946 {
02947
02948 pv[i].setTOF ( tof );
02949
02950 }
02951
02952 for(i=0; i<vecLen; i++)
02953 {
02954 if(pv[i].getName() == "KaonZero" || pv[i].getName() == "AntiKaonZero")
02955 {
02956 pvmx[0] = pv[i];
02957 if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
02958 else pv[i].setDefinition("KaonZeroLong");
02959 pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
02960 }
02961 }
02962
02963 successful = true;
02964 delete [] pvmx;
02965 delete [] tempV;
02966 return;
02967 }
02968
02969 void
02970 G4HEInelastic::MediumEnergyCascading(G4bool& successful,
02971 G4HEVector pv[],
02972 G4int& vecLen,
02973 G4double& excitationEnergyGNP,
02974 G4double& excitationEnergyDTA,
02975 const G4HEVector& incidentParticle,
02976 const G4HEVector& targetParticle,
02977 G4double atomicWeight,
02978 G4double atomicNumber)
02979 {
02980
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990 G4int protonCode = Proton.getCode();
02991 G4double protonMass = Proton.getMass();
02992 G4int neutronCode = Neutron.getCode();
02993 G4double kaonPlusMass = KaonPlus.getMass();
02994 G4int kaonPlusCode = KaonPlus.getCode();
02995 G4int kaonMinusCode = KaonMinus.getCode();
02996 G4int kaonZeroSCode = KaonZeroShort.getCode();
02997 G4int kaonZeroLCode = KaonZeroLong.getCode();
02998 G4int kaonZeroCode = KaonZero.getCode();
02999 G4int antiKaonZeroCode = AntiKaonZero.getCode();
03000 G4int pionPlusCode = PionPlus.getCode();
03001 G4int pionZeroCode = PionZero.getCode();
03002 G4int pionMinusCode = PionMinus.getCode();
03003 G4String mesonType = PionPlus.getType();
03004 G4String baryonType = Proton.getType();
03005 G4String antiBaryonType = AntiProton.getType();
03006
03007 G4double targetMass = targetParticle.getMass();
03008
03009 G4int incidentCode = incidentParticle.getCode();
03010 G4double incidentMass = incidentParticle.getMass();
03011 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
03012 G4double incidentEnergy = incidentParticle.getEnergy();
03013 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
03014 G4String incidentType = incidentParticle.getType();
03015
03016 G4double incidentTOF = 0.;
03017
03018
03019
03020 G4int i, j, l;
03021
03022 if(verboseLevel > 1)
03023 G4cout << " G4HEInelastic::MediumEnergyCascading " << G4endl;
03024
03025
03026
03027 G4bool annihilation = false;
03028 if (incidentCode < 0 && incidentType == antiBaryonType &&
03029 pv[0].getType() != antiBaryonType &&
03030 pv[1].getType() != antiBaryonType) {
03031 annihilation = true;
03032 }
03033
03034 successful = false;
03035
03036 G4double twsup[] = { 1., 1., 0.7, 0.5, 0.3, 0.2, 0.1, 0.0 };
03037
03038 if(annihilation) goto start;
03039 if(vecLen >= 8) goto start;
03040 if(incidentKineticEnergy < 1.) return;
03041 if( ( incidentCode == kaonPlusCode || incidentCode == kaonMinusCode
03042 || incidentCode == kaonZeroCode || incidentCode == antiKaonZeroCode
03043 || incidentCode == kaonZeroSCode || incidentCode == kaonZeroLCode )
03044 && ( G4UniformRand() < 0.5)) goto start;
03045 if(G4UniformRand() > twsup[vecLen-1]) goto start;
03046 return;
03047
03048 start:
03049
03050 if (annihilation)
03051 {
03052 G4double ekcor = Amax( 1., 1./incidentKineticEnergy);
03053 incidentKineticEnergy = 2*targetMass + incidentKineticEnergy*(1.+ekcor/atomicWeight);
03054 G4double excitation = NuclearExcitation(incidentKineticEnergy,
03055 atomicWeight,
03056 atomicNumber,
03057 excitationEnergyGNP,
03058 excitationEnergyDTA);
03059 incidentKineticEnergy -= excitation;
03060 if (incidentKineticEnergy < excitationEnergyDTA) incidentKineticEnergy = 0.;
03061 incidentEnergy = incidentKineticEnergy + incidentMass;
03062 incidentTotalMomentum =
03063 std::sqrt( Amax(0., incidentEnergy*incidentEnergy - incidentMass*incidentMass));
03064 }
03065
03066 G4HEVector pTemp;
03067 for(i = 2; i<vecLen; i++)
03068 {
03069 j = Imin(vecLen-1, (G4int)(2. + G4UniformRand()*(vecLen-2)));
03070 pTemp = pv[j];
03071 pv[j] = pv[i];
03072 pv[i] = pTemp;
03073 }
03074
03075
03076
03077
03078
03079 if ((incidentCode==kaonPlusCode || incidentCode==kaonMinusCode ||
03080 incidentCode==kaonZeroCode || incidentCode==antiKaonZeroCode ||
03081 incidentCode==kaonZeroSCode || incidentCode==kaonZeroLCode)
03082 && (G4UniformRand()>0.7) ) {
03083 pTemp = pv[1];
03084 pv[1] = pv[0];
03085 pv[0] = pTemp;
03086 }
03087
03088
03089
03090
03091
03092
03093
03094 G4int lead = 0;
03095 G4HEVector leadParticle;
03096 if ((incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
03097 && (incidentCode != neutronCode) ) {
03098 G4double pMass = pv[0].getMass();
03099 G4int pCode = pv[0].getCode();
03100 if ((pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
03101 && (pCode != neutronCode) ) {
03102 lead = pCode;
03103 leadParticle = pv[0];
03104 } else {
03105 pMass = pv[1].getMass();
03106 pCode = pv[1].getCode();
03107 if ((pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
03108 && (pCode != neutronCode) ) {
03109 lead = pCode;
03110 leadParticle = pv[1];
03111 }
03112 }
03113 }
03114
03115
03116
03117
03118 G4HEVector pvI = incidentParticle;
03119 pvI.setSide( 1 );
03120
03121 G4HEVector pvT = targetParticle;
03122 pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
03123 pvT.setSide( -1 );
03124 pvT.setTOF( -1.);
03125
03126 G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass)+sqr(targetMass)
03127 +2.0*targetMass*incidentEnergy );
03128
03129 G4double tavai1 = centerOfMassEnergy/2.0 - incidentMass;
03130 G4double tavai2 = centerOfMassEnergy/2.0 - targetMass;
03131
03132 G4int ntb = 1;
03133 for (i = 0; i < vecLen; i++) {
03134 if (i == 0) {
03135 pv[i].setSide(1);
03136 } else if (i == 1) {
03137 pv[i].setSide(-1);
03138 } else {
03139 if (G4UniformRand() < 0.5) {
03140 pv[i].setSide(-1);
03141 ntb++;
03142 } else pv[i].setSide(1);
03143 }
03144 pv[i].setTOF(incidentTOF);
03145 }
03146
03147 G4double tb = 2. * ntb;
03148 if (centerOfMassEnergy < (2. + G4UniformRand()))
03149 tb = (2. * ntb + vecLen)/2.;
03150
03151 if (verboseLevel > 1) {
03152 G4cout << " pv Vector after Randomization " << vecLen << G4endl;
03153 pvI.Print(-1);
03154 pvT.Print(-1);
03155 for (i=0; i < vecLen ; i++) pv[i].Print(i);
03156 }
03157
03158
03159
03160
03161
03162
03163 G4double ss, xtarg, ran;
03164 ss = centerOfMassEnergy*centerOfMassEnergy;
03165 xtarg = Amax( 0.01, Amin( 0.75, 0.312 + 0.200 * std::log(std::log(ss))
03166 + std::pow(ss,1.5)/6000.0 )
03167 *(std::pow(atomicWeight, 0.33) - 1.0) * tb);
03168
03169 G4int ntarg = Poisson(xtarg);
03170 G4int targ = 0;
03171
03172 if (ntarg > 0) {
03173 G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
03174 G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
03175 G4int momentumBin = 0;
03176 while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin]) )
03177 momentumBin++;
03178 momentumBin = Imin( 5, momentumBin );
03179
03180
03181
03182
03183 for( i=0; i<ntarg; i++ )
03184 {
03185 if( G4UniformRand() < nucsup[momentumBin] )
03186 {
03187 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
03188 pv[vecLen].setDefinition( "Proton" );
03189 else
03190 pv[vecLen].setDefinition( "Neutron" );
03191 targ++;
03192 }
03193 else
03194 {
03195 ran = G4UniformRand();
03196 if( ran < 0.33333 )
03197 pv[vecLen].setDefinition( "PionPlus");
03198 else if( ran < 0.66667 )
03199 pv[vecLen].setDefinition( "PionZero");
03200 else
03201 pv[vecLen].setDefinition( "PionMinus" );
03202 }
03203 pv[vecLen].setSide( -2 );
03204 pv[vecLen].setFlag( true );
03205 pv[vecLen].setTOF( incidentTOF );
03206 vecLen++;
03207 }
03208 }
03209
03210
03211
03212
03213 G4int is, iskip;
03214 tavai1 = centerOfMassEnergy/2.;
03215 G4int iavai1 = 0;
03216
03217 for (i = 0; i < vecLen; i++)
03218 {
03219 if (pv[i].getSide() > 0)
03220 {
03221 tavai1 -= pv[i].getMass();
03222 iavai1++;
03223 }
03224 }
03225 if ( iavai1 == 0) return;
03226
03227 while( tavai1 <= 0.0 )
03228 {
03229 iskip = G4int(G4UniformRand()*iavai1) + 1;
03230 is = 0;
03231 for( i=vecLen-1; i>=0; i-- )
03232 {
03233 if( pv[i].getSide() > 0 )
03234 {
03235 if (++is == iskip)
03236 {
03237 tavai1 += pv[i].getMass();
03238 iavai1--;
03239 if ( i != vecLen-1)
03240 {
03241 for( j=i; j<vecLen; j++ )
03242 {
03243 pv[j] = pv[j+1];
03244 }
03245 }
03246 if( --vecLen == 0 ) return;
03247 break;
03248 }
03249 }
03250 }
03251 }
03252
03253
03254 tavai2 = (targ+1)*centerOfMassEnergy/2.;
03255 G4int iavai2 = 0;
03256
03257 for (i = 0; i < vecLen; i++)
03258 {
03259 if (pv[i].getSide() < 0)
03260 {
03261 tavai2 -= pv[i].getMass();
03262 iavai2++;
03263 }
03264 }
03265 if (iavai2 == 0) return;
03266
03267 while( tavai2 <= 0.0 )
03268 {
03269 iskip = G4int(G4UniformRand()*iavai2) + 1;
03270 is = 0;
03271 for( i = vecLen-1; i >= 0; i-- )
03272 {
03273 if( pv[i].getSide() < 0 )
03274 {
03275 if( ++is == iskip )
03276 {
03277 tavai2 += pv[i].getMass();
03278 iavai2--;
03279 if (pv[i].getSide() == -2) ntarg--;
03280 if (i != vecLen-1)
03281 {
03282 for( j=i; j<vecLen; j++)
03283 {
03284 pv[j] = pv[j+1];
03285 }
03286 }
03287 if (--vecLen == 0) return;
03288 break;
03289 }
03290 }
03291 }
03292 }
03293
03294 if (verboseLevel > 1) {
03295 G4cout << " pv Vector after Energy checks " << vecLen << " "
03296 << tavai1 << " " << iavai1 << " " << tavai2 << " "
03297 << iavai2 << " " << ntarg << G4endl;
03298 pvI.Print(-1);
03299 pvT.Print(-1);
03300 for (i=0; i < vecLen ; i++) pv[i].Print(i);
03301 }
03302
03303
03304
03305 G4HEVector* pvmx = new G4HEVector [10];
03306
03307 pvmx[0].setMass( incidentMass );
03308 pvmx[0].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
03309 pvmx[1].setMass( protonMass);
03310 pvmx[1].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
03311 pvmx[3].setMass( protonMass*(1+targ));
03312 pvmx[3].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
03313 pvmx[4].setZero();
03314 pvmx[5].setZero();
03315 pvmx[7].setZero();
03316 pvmx[8].setZero();
03317 pvmx[8].setMomentum( 1.0, 0.0 );
03318 pvmx[2].Add( pvmx[0], pvmx[1] );
03319 pvmx[3].Add( pvmx[3], pvmx[0] );
03320 pvmx[0].Lor( pvmx[0], pvmx[2] );
03321 pvmx[1].Lor( pvmx[1], pvmx[2] );
03322
03323 if (verboseLevel > 1) {
03324 G4cout << " General Vectors after Definition " << G4endl;
03325 for (i=0; i<10; i++) pvmx[i].Print(i);
03326 }
03327
03328
03329
03330
03331
03332 G4double dndl[20];
03333 G4double binl[20];
03334 G4double pvMass, pvEnergy;
03335 G4int pvCode;
03336 G4double aspar, pt, phi, et, xval;
03337 G4double ekin = 0.;
03338 G4double ekin1 = 0.;
03339 G4double ekin2 = 0.;
03340 phi = G4UniformRand()*twopi;
03341 G4int npg = 0;
03342 G4int targ1 = 0;
03343 for( i=vecLen-1; i>=0; i-- )
03344 {
03345 if( (pv[i].getSide() == -2) || (i == 1) )
03346 {
03347 if ( pv[i].getType() == baryonType ||
03348 pv[i].getType() == antiBaryonType)
03349 {
03350 if( ++npg < 19 )
03351 {
03352 pv[i].setSide( -3 );
03353 targ1++;
03354 continue;
03355 }
03356 }
03357 }
03358
03359
03360
03361
03362
03363 G4double maspar[] = { 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.75, 0.20};
03364 G4double bp[] = { 3.50, 3.50, 3.50, 6.00, 5.00, 4.00, 3.50, 3.50};
03365 G4double ptex[] = { 1.70, 1.70, 1.50, 1.70, 1.40, 1.20, 1.70, 1.20};
03366
03367
03368
03369
03370 pvMass = pv[i].getMass();
03371 j = 2;
03372 if ( pv[i].getType() == mesonType ) j = 1;
03373 if ( pv[i].getMass() < 0.4 ) j = 0;
03374 if ( i <= 1 ) j += 3;
03375 if (pv[i].getSide() <= -2) j = 6;
03376 if (j == 6 && (pv[i].getType() == baryonType || pv[i].getType()==antiBaryonType) ) j = 7;
03377 pt = Amax(0.001, std::sqrt(std::pow(-std::log(1.-G4UniformRand())/bp[j],ptex[j])));
03378 aspar = maspar[j];
03379 phi = G4UniformRand()*twopi;
03380 pv[i].setMomentum( pt*std::cos(phi), pt*std::sin(phi) );
03381
03382 for( j=0; j<20; j++ ) binl[j] = j/(19.*pt);
03383
03384 if( pv[i].getSide() > 0 )
03385 et = pvmx[0].getEnergy();
03386 else
03387 et = pvmx[1].getEnergy();
03388
03389 dndl[0] = 0.0;
03390
03391
03392
03393 G4int outerCounter = 0, innerCounter = 0;
03394 G4bool eliminateThisParticle = true;
03395 G4bool resetEnergies = true;
03396 while( ++outerCounter < 3 )
03397 {
03398 for( l=1; l<20; l++ )
03399 {
03400 xval = (binl[l]+binl[l-1])/2.;
03401 if( xval > 1./pt )
03402 dndl[l] = dndl[l-1];
03403 else
03404 dndl[l] = dndl[l-1] +
03405 aspar/std::sqrt( std::pow((1.+aspar*xval*aspar*xval),3) ) *
03406 (binl[l]-binl[l-1]) * et /
03407 std::sqrt( pt*xval*et*pt*xval*et + pt*pt + pvMass*pvMass );
03408 }
03409
03410
03411
03412 innerCounter = 0;
03413 while( ++innerCounter < 7 )
03414 {
03415 l = 1;
03416 ran = G4UniformRand()*dndl[19];
03417 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
03418 l = Imin( 19, l );
03419 xval = Amin( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
03420 if( pv[i].getSide() < 0 ) xval *= -1.;
03421 pv[i].setMomentumAndUpdate( xval*et );
03422 pvEnergy = pv[i].getEnergy();
03423 if( pv[i].getSide() > 0 )
03424 {
03425 if ( i < 2 )
03426 {
03427 ekin = tavai1 - ekin1;
03428 if (ekin < 0.) ekin = 0.04*std::fabs(normal());
03429 G4double pp1 = pv[i].Length();
03430 if (pp1 >= 1.e-6)
03431 {
03432 G4double pp = std::sqrt(ekin*(ekin+2*pvMass));
03433 pp = Amax(0.,pp*pp-pt*pt);
03434 pp = std::sqrt(pp)*pv[i].getSide()/std::fabs(G4double(pv[i].getSide()));
03435 pv[i].setMomentumAndUpdate( pp );
03436 }
03437 else
03438 {
03439 pv[i].setMomentum(0.,0.,0.);
03440 pv[i].setKineticEnergyAndUpdate( ekin);
03441 }
03442 pvmx[4].Add( pvmx[4], pv[i]);
03443 outerCounter = 2;
03444 resetEnergies = false;
03445 eliminateThisParticle = false;
03446 break;
03447 }
03448 else if( (ekin1+pvEnergy-pvMass) < 0.95*tavai1 )
03449 {
03450 pvmx[4].Add( pvmx[4], pv[i] );
03451 ekin1 += pvEnergy - pvMass;
03452 pvmx[6].Add( pvmx[4], pvmx[5] );
03453 pvmx[6].setMomentum( 0.0 );
03454 outerCounter = 2;
03455 eliminateThisParticle = false;
03456 resetEnergies = false;
03457 break;
03458 }
03459 if( innerCounter > 5 ) break;
03460
03461 if( tavai2 >= pvMass )
03462 {
03463 pv[i].setSide( -1 );
03464 tavai1 += pvMass;
03465 tavai2 -= pvMass;
03466 iavai2++;
03467 }
03468 }
03469 else
03470 {
03471 xval = Amin(0.999,0.95+0.05*targ/20.0);
03472 if( (ekin2+pvEnergy-pvMass) < xval*tavai2 )
03473 {
03474 pvmx[5].Add( pvmx[5], pv[i] );
03475 ekin2 += pvEnergy - pvMass;
03476 pvmx[6].Add( pvmx[4], pvmx[5] );
03477 pvmx[6].setMomentum( 0.0 );
03478 outerCounter = 2;
03479 eliminateThisParticle = false;
03480 resetEnergies = false;
03481 break;
03482 }
03483 if( innerCounter > 5 )break;
03484
03485 if( tavai1 >= pvMass )
03486 {
03487 pv[i].setSide( 1 );
03488 tavai1 -= pvMass;
03489 tavai2 += pvMass;
03490 iavai2--;
03491 }
03492 }
03493 pv[i].setMomentum( pv[i].getMomentum().x() * 0.9,
03494 pv[i].getMomentum().y() * 0.9);
03495 pt *= 0.9;
03496 dndl[19] *= 0.9;
03497 }
03498
03499 if (resetEnergies)
03500 {
03501 ekin1 = 0.0;
03502 ekin2 = 0.0;
03503 pvmx[4].setZero();
03504 pvmx[5].setZero();
03505 if (verboseLevel > 1)
03506 G4cout << " Reset energies for index " << i << G4endl;
03507 for( l=i+1; l < vecLen; l++ )
03508 {
03509 if( (pv[l].getMass() < protonMass) || (pv[l].getSide() > 0) )
03510 {
03511 pvEnergy = Amax( pv[l].getMass(), 0.95*pv[l].getEnergy()
03512 + 0.05*pv[l].getMass() );
03513 pv[l].setEnergyAndUpdate( pvEnergy );
03514 if( pv[l].getSide() > 0)
03515 {
03516 ekin1 += pv[l].getKineticEnergy();
03517 pvmx[4].Add( pvmx[4], pv[l] );
03518 }
03519 else
03520 {
03521 ekin2 += pv[l].getKineticEnergy();
03522 pvmx[5].Add( pvmx[5], pv[l] );
03523 }
03524 }
03525 }
03526 }
03527 }
03528
03529 if( eliminateThisParticle )
03530 {
03531 if (verboseLevel > 1)
03532 {
03533 G4cout << " Eliminate particle with index " << i << G4endl;
03534 pv[i].Print(i);
03535 }
03536 for( j=i; j < vecLen; j++ )
03537 {
03538 pv[j] = pv[j+1];
03539 }
03540 vecLen--;
03541 if (vecLen < 2) {
03542 delete [] pvmx;
03543 return;
03544 }
03545 i++;
03546 pvmx[6].Add( pvmx[4], pvmx[5] );
03547 pvmx[6].setMomentum( 0.0 );
03548 }
03549 }
03550 if (verboseLevel > 1)
03551 { G4cout << " pv Vector after lambda fragmentation " << vecLen << G4endl;
03552 pvI.Print(-1);
03553 pvT.Print(-1);
03554 for (i=0; i < vecLen ; i++) pv[i].Print(i);
03555 for (i=0; i < 10; i++) pvmx[i].Print(i);
03556 }
03557
03558
03559
03560 pvmx[6].Lor( pvmx[3], pvmx[2] );
03561 pvmx[6].Sub( pvmx[6], pvmx[4] );
03562 pvmx[6].Sub( pvmx[6], pvmx[5] );
03563 if (verboseLevel > 1) pvmx[6].Print(6);
03564
03565 npg = 0;
03566 G4double rmb0 = 0.;
03567 G4double rmb;
03568 G4double wgt;
03569 G4bool constantCrossSection = true;
03570 for (i = 0; i < vecLen; i++)
03571 {
03572 if(pv[i].getSide() == -3)
03573 {
03574 npg++;
03575 rmb0 += pv[i].getMass();
03576 }
03577 }
03578 if( targ1 == 1 || npg < 2)
03579 {
03580 ekin = Amin( tavai2-ekin2, centerOfMassEnergy/2.0-protonMass );
03581 if( ekin < 0.04 ) ekin = 0.04 * std::fabs( normal() );
03582 G4double pp = std::sqrt(ekin*(ekin+2*pv[1].getMass()));
03583 G4double pp1 = pvmx[6].Length();
03584 if(pp1 < 1.e-6)
03585 {
03586 pv[1].setKineticEnergyAndUpdate(ekin);
03587 }
03588 else
03589 {
03590 pv[1].setMomentum(pvmx[6].getMomentum());
03591 pv[1].SmulAndUpdate(pv[1],pp/pp1);
03592 }
03593 pvmx[5].Add( pvmx[5], pv[1] );
03594 }
03595 else
03596 {
03597 G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
03598 G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
03599
03600 G4int tempCount = Imin( 5, targ1 ) - 1;
03601
03602 rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()), cpar[tempCount])/gpar[tempCount];
03603 pvEnergy = pvmx[6].getEnergy();
03604 if ( rmb > pvEnergy ) rmb = pvEnergy;
03605 pvmx[6].setMass( rmb );
03606 pvmx[6].setEnergyAndUpdate( pvEnergy );
03607 pvmx[6].Smul( pvmx[6], -1. );
03608 if (verboseLevel > 1) {
03609 G4cout << " General Vectors before input to NBodyPhaseSpace "
03610 << targ1 << " " << tempCount << " " << rmb0 << " "
03611 << rmb << " " << pvEnergy << G4endl;
03612 for (i=0; i<10; i++) pvmx[i].Print(i);
03613 }
03614
03615
03616
03617 G4HEVector* tempV = new G4HEVector[18];
03618 npg = 0;
03619 for( i=0; i < vecLen; i++ )
03620 {
03621 if( pv[i].getSide() == -3 ) tempV[npg++] = pv[i];
03622 }
03623
03624 wgt = NBodyPhaseSpace( pvmx[6].getMass(), constantCrossSection, tempV, npg );
03625
03626 npg = 0;
03627 for( i=0; i < vecLen; i++ )
03628 {
03629 if( pv[i].getSide() == -3 )
03630 {
03631 pv[i].setMomentum( tempV[npg++].getMomentum());
03632 pv[i].SmulAndUpdate( pv[i], 1.);
03633 pv[i].Lor( pv[i], pvmx[6] );
03634 pvmx[5].Add( pvmx[5], pv[i] );
03635 }
03636 }
03637 delete [] tempV;
03638 }
03639 if( vecLen <= 2 )
03640 {
03641 successful = false;
03642 return;
03643 }
03644
03645
03646
03647 targ = 0;
03648 for( i=0; i < vecLen; i++ )
03649 {
03650 if( pv[i].getType() == baryonType )targ++;
03651 if( pv[i].getType() == antiBaryonType )targ++;
03652 pv[i].Lor( pv[i], pvmx[1] );
03653 }
03654 targ = Imax( 1, targ );
03655
03656 G4bool dum(0);
03657 if( lead )
03658 {
03659 for( i=0; i<vecLen; i++ )
03660 {
03661 if( pv[i].getCode() == lead )
03662 {
03663 dum = false;
03664 break;
03665 }
03666 }
03667 if( dum )
03668 {
03669 i = 0;
03670
03671 if( ( (leadParticle.getType() == baryonType ||
03672 leadParticle.getType() == antiBaryonType)
03673 && (pv[1].getType() == baryonType ||
03674 pv[1].getType() == antiBaryonType))
03675 || ( (leadParticle.getType() == mesonType)
03676 && (pv[1].getType() == mesonType)))
03677 {
03678 i = 1;
03679 }
03680 ekin = pv[i].getKineticEnergy();
03681 pv[i] = leadParticle;
03682 if( pv[i].getFlag() )
03683 pv[i].setTOF( -1.0 );
03684 else
03685 pv[i].setTOF( 1.0 );
03686 pv[i].setKineticEnergyAndUpdate( ekin );
03687 }
03688 }
03689
03690 pvmx[3].setMass( incidentMass);
03691 pvmx[3].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
03692
03693 G4double ekin0 = pvmx[3].getKineticEnergy();
03694
03695 pvmx[4].setMass ( protonMass * targ);
03696 pvmx[4].setEnergy( protonMass * targ);
03697 pvmx[4].setMomentum(0.,0.,0.);
03698 pvmx[4].setKineticEnergy(0.);
03699
03700 ekin = pvmx[3].getEnergy() + pvmx[4].getEnergy();
03701
03702 pvmx[5].Add( pvmx[3], pvmx[4] );
03703 pvmx[3].Lor( pvmx[3], pvmx[5] );
03704 pvmx[4].Lor( pvmx[4], pvmx[5] );
03705
03706 G4double tecm = pvmx[3].getEnergy() + pvmx[4].getEnergy();
03707
03708 pvmx[7].setZero();
03709
03710 ekin1 = 0.0;
03711 G4double teta;
03712
03713 for( i=0; i < vecLen; i++ )
03714 {
03715 pvmx[7].Add( pvmx[7], pv[i] );
03716 ekin1 += pv[i].getKineticEnergy();
03717 ekin -= pv[i].getMass();
03718 }
03719
03720 if( vecLen > 1 && vecLen < 19 )
03721 {
03722 constantCrossSection = true;
03723 G4HEVector pw[18];
03724 for(i=0;i<vecLen;i++) pw[i] = pv[i];
03725 wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
03726 ekin = 0.0;
03727 for( i=0; i < vecLen; i++ )
03728 {
03729 pvmx[6].setMass( pw[i].getMass());
03730 pvmx[6].setMomentum( pw[i].getMomentum() );
03731 pvmx[6].SmulAndUpdate( pvmx[6], 1.);
03732 pvmx[6].Lor( pvmx[6], pvmx[4] );
03733 ekin += pvmx[6].getKineticEnergy();
03734 }
03735 teta = pvmx[7].Ang( pvmx[3] );
03736 if (verboseLevel > 1)
03737 G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
03738 << " " << ekin1 << " " << ekin << G4endl;
03739 }
03740
03741 if( ekin1 != 0.0 )
03742 {
03743 pvmx[6].setZero();
03744 wgt = ekin/ekin1;
03745 ekin1 = 0.;
03746 for( i=0; i < vecLen; i++ )
03747 {
03748 pvMass = pv[i].getMass();
03749 ekin = pv[i].getKineticEnergy() * wgt;
03750 pv[i].setKineticEnergyAndUpdate( ekin );
03751 ekin1 += ekin;
03752 pvmx[6].Add( pvmx[6], pv[i] );
03753 }
03754 teta = pvmx[6].Ang( pvmx[3] );
03755 if (verboseLevel > 1)
03756 G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
03757 << ekin1 << G4endl;
03758 }
03759
03760
03761
03762 G4double ry = G4UniformRand();
03763 G4double rz = G4UniformRand();
03764 G4double rx = twopi*rz;
03765 G4double a1 = std::sqrt(-2.0*std::log(ry));
03766 G4double rantarg1 = a1*std::cos(rx)*0.02*targ/G4double(vecLen);
03767 G4double rantarg2 = a1*std::sin(rx)*0.02*targ/G4double(vecLen);
03768
03769 for (i = 0; i < vecLen; i++)
03770 pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
03771 pv[i].getMomentum().y()+rantarg2 );
03772
03773 if (verboseLevel > 1) {
03774 pvmx[6].setZero();
03775 for (i = 0; i < vecLen; i++) pvmx[6].Add( pvmx[6], pv[i] );
03776 teta = pvmx[6].Ang( pvmx[3] );
03777 G4cout << " After smearing " << teta << G4endl;
03778 }
03779
03780
03781
03782
03783
03784
03785
03786
03787 G4double dekin = 0.0;
03788 G4int npions = 0;
03789 G4double ek1 = 0.0;
03790 G4double alekw, xxh;
03791 G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
03792 G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.00};
03793 G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
03794
03795 for (i = 0; i < vecLen; i++) {
03796 pv[i].Defs1( pv[i], pvI );
03797 if (atomicWeight > 1.5) {
03798 ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa*( 1. + 0.5*normal()));
03799 alekw = std::log( incidentKineticEnergy );
03800 xxh = 1.;
03801 if (incidentCode == pionPlusCode || incidentCode == pionMinusCode) {
03802 if (pv[i].getCode() == pionZeroCode) {
03803 if (G4UniformRand() < std::log(atomicWeight)) {
03804 if (alekw > alem[0]) {
03805 G4int jmax = 1;
03806 for (j = 1; j < 8; j++) {
03807 if (alekw < alem[j]) {
03808 jmax = j;
03809 break;
03810 }
03811 }
03812 xxh = (val0[jmax] - val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alekw
03813 + val0[jmax-1] - (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alem[jmax-1];
03814 xxh = 1. - xxh;
03815 }
03816 }
03817 }
03818 }
03819 dekin += ekin*(1.-xxh);
03820 ekin *= xxh;
03821 pv[i].setKineticEnergyAndUpdate( ekin );
03822 pvCode = pv[i].getCode();
03823 if ((pvCode == pionPlusCode) ||
03824 (pvCode == pionMinusCode) ||
03825 (pvCode == pionZeroCode)) {
03826 npions += 1;
03827 ek1 += ekin;
03828 }
03829 }
03830 }
03831
03832 if( (ek1 > 0.0) && (npions > 0) )
03833 {
03834 dekin = 1.+dekin/ek1;
03835 for (i = 0; i < vecLen; i++)
03836 {
03837 pvCode = pv[i].getCode();
03838 if((pvCode == pionPlusCode) || (pvCode == pionMinusCode) || (pvCode == pionZeroCode))
03839 {
03840 ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
03841 pv[i].setKineticEnergyAndUpdate( ekin );
03842 }
03843 }
03844 }
03845 if (verboseLevel > 1)
03846 { G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
03847 for (i=0; i<vecLen; i++) pv[i].Print(i);
03848 }
03849
03850
03851
03852
03853
03854 if (verboseLevel > 1) G4cout << " Evaporation " << atomicWeight << " " <<
03855 excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
03856
03857 if( atomicWeight > 1.5 )
03858 {
03859
03860 G4double sprob, cost, sint, pp, eka;
03861 G4int spall(0), nbl(0);
03862
03863
03864 if( incidentKineticEnergy < 5.0 )
03865 sprob = 0.0;
03866 else
03867
03868 sprob = Amin(1., 0.000314*atomicWeight*std::log(incidentKineticEnergy-4.));
03869
03870
03871
03872 if( excitationEnergyGNP >= 0.001 )
03873 {
03874
03875
03876
03877 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyGNP/
03878 (excitationEnergyGNP+excitationEnergyDTA));
03879 if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
03880 if (verboseLevel > 1)
03881 G4cout << " evaporation " << targ << " " << nbl << " "
03882 << sprob << G4endl;
03883 spall = targ;
03884 if( nbl > 0)
03885 {
03886 ekin = excitationEnergyGNP/nbl;
03887 ekin2 = 0.0;
03888 for( i=0; i<nbl; i++ )
03889 {
03890 if( G4UniformRand() < sprob ) continue;
03891 if( ekin2 > excitationEnergyGNP) break;
03892 ran = G4UniformRand();
03893 ekin1 = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
03894 if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
03895 ekin2 += ekin1;
03896 if( ekin2 > excitationEnergyGNP)
03897 ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
03898 if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
03899 pv[vecLen].setDefinition( "Proton");
03900 else
03901 pv[vecLen].setDefinition( "Neutron");
03902 spall++;
03903 cost = G4UniformRand() * 2.0 - 1.0;
03904 sint = std::sqrt(std::fabs(1.0-cost*cost));
03905 phi = twopi * G4UniformRand();
03906 pv[vecLen].setFlag( true );
03907 pv[vecLen].setSide( -4 );
03908 pvMass = pv[vecLen].getMass();
03909 pv[vecLen].setTOF( 1.0 );
03910 pvEnergy = ekin1 + pvMass;
03911 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
03912 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
03913 pp*sint*std::cos(phi),
03914 pp*cost );
03915 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
03916 vecLen++;
03917 }
03918 if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
03919 {
03920 G4int ika, kk = 0;
03921 eka = incidentKineticEnergy;
03922 if( eka > 1.0 )eka *= eka;
03923 eka = Amax( 0.1, eka );
03924 ika = G4int(3.6*std::exp((atomicNumber*atomicNumber
03925 /atomicWeight-35.56)/6.45)/eka);
03926 if( ika > 0 )
03927 {
03928 for( i=(vecLen-1); i>=0; i-- )
03929 {
03930 if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
03931 {
03932 pTemp = pv[i];
03933 pv[i].setDefinition( "Neutron");
03934 pv[i].setMomentumAndUpdate(pTemp.getMomentum());
03935 if (verboseLevel > 1) pv[i].Print(i);
03936 if( ++kk > ika ) break;
03937 }
03938 }
03939 }
03940 }
03941 }
03942 }
03943
03944
03945
03946
03947 if( excitationEnergyDTA >= 0.001 )
03948 {
03949 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyDTA
03950 /(excitationEnergyGNP+excitationEnergyDTA));
03951
03952
03953
03954 if( nbl > 0 )
03955 {
03956 ekin = excitationEnergyDTA/nbl;
03957 ekin2 = 0.0;
03958 for( i=0; i<nbl; i++ )
03959 {
03960 if( G4UniformRand() < sprob ) continue;
03961 if( ekin2 > excitationEnergyDTA) break;
03962 ran = G4UniformRand();
03963 ekin1 = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
03964 if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
03965 ekin2 += ekin1;
03966 if( ekin2 > excitationEnergyDTA)
03967 ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
03968 cost = G4UniformRand()*2.0 - 1.0;
03969 sint = std::sqrt(std::fabs(1.0-cost*cost));
03970 phi = twopi*G4UniformRand();
03971 ran = G4UniformRand();
03972 if( ran <= 0.60 )
03973 pv[vecLen].setDefinition( "Deuteron");
03974 else if (ran <= 0.90)
03975 pv[vecLen].setDefinition( "Triton");
03976 else
03977 pv[vecLen].setDefinition( "Alpha");
03978 spall += (int)(pv[vecLen].getMass() * 1.066);
03979 if( spall > atomicWeight ) break;
03980 pv[vecLen].setFlag( true );
03981 pv[vecLen].setSide( -4 );
03982 pvMass = pv[vecLen].getMass();
03983 pv[vecLen].setSide( pv[vecLen].getCode());
03984 pv[vecLen].setTOF( 1.0 );
03985 pvEnergy = pvMass + ekin1;
03986 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
03987 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
03988 pp*sint*std::cos(phi),
03989 pp*cost );
03990 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
03991 vecLen++;
03992 }
03993 }
03994 }
03995 }
03996 if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
03997 {
03998 for( i=0; i<vecLen; i++ )
03999 {
04000 G4double etb = pv[i].getKineticEnergy();
04001 if( etb >= incidentKineticEnergy )
04002 pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
04003 }
04004 }
04005
04006
04007
04008 G4double tof = incidentTOF;
04009 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
04010 && (incidentKineticEnergy <= 0.2) )
04011 tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
04012 for ( i=0; i < vecLen; i++)
04013 {
04014
04015 pv[i].setTOF ( tof );
04016
04017 }
04018
04019 for(i=0; i<vecLen; i++)
04020 {
04021 if(pv[i].getName() == "KaonZero" || pv[i].getName() == "AntiKaonZero")
04022 {
04023 pvmx[0] = pv[i];
04024 if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
04025 else pv[i].setDefinition("KaonZeroLong");
04026 pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
04027 }
04028 }
04029
04030 successful = true;
04031 delete [] pvmx;
04032 return;
04033 }
04034
04035 void
04036 G4HEInelastic::MediumEnergyClusterProduction(G4bool& successful,
04037 G4HEVector pv[],
04038 G4int& vecLen,
04039 G4double& excitationEnergyGNP,
04040 G4double& excitationEnergyDTA,
04041 const G4HEVector& incidentParticle,
04042 const G4HEVector& targetParticle,
04043 G4double atomicWeight,
04044 G4double atomicNumber)
04045 {
04046
04047
04048
04049
04050
04051
04052
04053 G4int protonCode = Proton.getCode();
04054 G4double protonMass = Proton.getMass();
04055 G4int neutronCode = Neutron.getCode();
04056 G4double kaonPlusMass = KaonPlus.getMass();
04057 G4int pionPlusCode = PionPlus.getCode();
04058 G4int pionZeroCode = PionZero.getCode();
04059 G4int pionMinusCode = PionMinus.getCode();
04060 G4String mesonType = PionPlus.getType();
04061 G4String baryonType = Proton.getType();
04062 G4String antiBaryonType = AntiProton.getType();
04063
04064 G4double targetMass = targetParticle.getMass();
04065
04066 G4int incidentCode = incidentParticle.getCode();
04067 G4double incidentMass = incidentParticle.getMass();
04068 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
04069 G4double incidentEnergy = incidentParticle.getEnergy();
04070 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
04071 G4String incidentType = incidentParticle.getType();
04072
04073 G4double incidentTOF = 0.;
04074
04075
04076
04077 G4int i, j;
04078
04079 if (verboseLevel > 1)
04080 G4cout << " G4HEInelastic::MediumEnergyClusterProduction " << G4endl;
04081
04082 if (incidentTotalMomentum < 0.01) {
04083 successful = false;
04084 return;
04085 }
04086 G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
04087 +2.*targetMass*incidentEnergy);
04088
04089 G4HEVector pvI = incidentParticle;
04090 pvI.setSide( 1 );
04091
04092 G4HEVector pvT = targetParticle;
04093 pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
04094 pvT.setSide( -1 );
04095 pvT.setTOF( -1.);
04096
04097
04098
04099
04100
04101 G4int targ = 0;
04102 G4int ifor = 0;
04103 G4int iback = 0;
04104 G4int pvCode;
04105 G4double pvMass, pvEnergy;
04106
04107 pv[0].setSide( 1 );
04108 pv[1].setSide( -1 );
04109 for(i = 0; i < vecLen; i++)
04110 {
04111 if (i > 1)
04112 {
04113 if( G4UniformRand() < 0.5)
04114 {
04115 pv[i].setSide( 1 );
04116 if (++ifor > 18)
04117 {
04118 pv[i].setSide( -1 );
04119 ifor--;
04120 iback++;
04121 }
04122 }
04123 else
04124 {
04125 pv[i].setSide( -1 );
04126 if (++iback > 18)
04127 {
04128 pv[i].setSide( 1 );
04129 ifor++;
04130 iback--;
04131 }
04132 }
04133 }
04134
04135 pvCode = pv[i].getCode();
04136
04137 if ( ( (incidentCode == protonCode) || (incidentCode == neutronCode)
04138 || (incidentType == mesonType) )
04139 && ( (pvCode == pionPlusCode) || (pvCode == pionMinusCode) )
04140 && ( (G4UniformRand() < (10.-incidentTotalMomentum)/6.) )
04141 && ( (G4UniformRand() < atomicWeight/300.) ) )
04142 {
04143 if (G4UniformRand() > atomicNumber/atomicWeight)
04144 pv[i].setDefinition( "Neutron");
04145 else
04146 pv[i].setDefinition( "Proton");
04147 targ++;
04148 }
04149 pv[i].setTOF( incidentTOF );
04150 }
04151 G4double tb = 2. * iback;
04152 if (centerOfMassEnergy < (2+G4UniformRand())) tb = (2.*iback + vecLen)/2.;
04153
04154 G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4};
04155
04156 G4double xtarg = Amax(0.01, (0.312+0.2*std::log(std::log(centerOfMassEnergy*centerOfMassEnergy)))
04157 * (std::pow(atomicWeight,0.33)-1.) * tb);
04158 G4int ntarg = Poisson(xtarg);
04159 if (ntarg > 0)
04160 {
04161 G4int ipx = Imin(4, (G4int)(incidentTotalMomentum/3.));
04162 for (i=0; i < ntarg; i++)
04163 {
04164 if (G4UniformRand() < nucsup[ipx] )
04165 {
04166 if (G4UniformRand() < (1.- atomicNumber/atomicWeight))
04167 pv[vecLen].setDefinition( "Neutron");
04168 else
04169 pv[vecLen].setDefinition( "Proton");
04170 targ++;
04171 }
04172 else
04173 {
04174 G4double ran = G4UniformRand();
04175 if (ran < 0.3333 )
04176 pv[vecLen].setDefinition( "PionPlus");
04177 else if (ran < 0.6666)
04178 pv[vecLen].setDefinition( "PionZero");
04179 else
04180 pv[vecLen].setDefinition( "PionMinus");
04181 }
04182 pv[vecLen].setSide( -2 );
04183 pv[vecLen].setFlag( true );
04184 pv[vecLen].setTOF( incidentTOF );
04185 vecLen++;
04186 }
04187 }
04188
04189
04190
04191
04192
04193
04194 G4int lead = 0;
04195 G4HEVector leadParticle;
04196 if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
04197 && (incidentCode != neutronCode) )
04198 {
04199 G4double pMass = pv[0].getMass();
04200 G4int pCode = pv[0].getCode();
04201 if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
04202 && (pCode != neutronCode) )
04203 {
04204 lead = pCode;
04205 leadParticle = pv[0];
04206 }
04207 else
04208 {
04209 pMass = pv[1].getMass();
04210 pCode = pv[1].getCode();
04211 if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
04212 && (pCode != neutronCode) )
04213 {
04214 lead = pCode;
04215 leadParticle = pv[1];
04216 }
04217 }
04218 }
04219
04220 if (verboseLevel > 1) {
04221 G4cout << " pv Vector after initialization " << vecLen << G4endl;
04222 pvI.Print(-1);
04223 pvT.Print(-1);
04224 for (i=0; i < vecLen ; i++) pv[i].Print(i);
04225 }
04226
04227 G4double tavai = 0.;
04228 for(i=0;i<vecLen;i++) if(pv[i].getSide() != -2) tavai += pv[i].getMass();
04229
04230 while (tavai > centerOfMassEnergy)
04231 {
04232 for (i=vecLen-1; i >= 0; i--)
04233 {
04234 if (pv[i].getSide() != -2)
04235 {
04236 tavai -= pv[i].getMass();
04237 if( i != vecLen-1)
04238 {
04239 for (j=i; j < vecLen; j++)
04240 {
04241 pv[j] = pv[j+1];
04242 }
04243 }
04244 if ( --vecLen < 2)
04245 {
04246 successful = false;
04247 return;
04248 }
04249 break;
04250 }
04251 }
04252 }
04253
04254
04255
04256
04257
04258
04259 G4double rmc0 = 0., rmd0 = 0., rme0 = 0.;
04260 G4int ntc = 0, ntd = 0, nte = 0;
04261
04262 for (i=0; i < vecLen; i++)
04263 {
04264 if(pv[i].getSide() > 0)
04265 {
04266 if(ntc < 17)
04267 {
04268 rmc0 += pv[i].getMass();
04269 ntc++;
04270 }
04271 else
04272 {
04273 if(ntd < 17)
04274 {
04275 pv[i].setSide(-1);
04276 rmd0 += pv[i].getMass();
04277 ntd++;
04278 }
04279 else
04280 {
04281 pv[i].setSide(-2);
04282 rme0 += pv[i].getMass();
04283 nte++;
04284 }
04285 }
04286 }
04287 else if (pv[i].getSide() == -1)
04288 {
04289 if(ntd < 17)
04290 {
04291 rmd0 += pv[i].getMass();
04292 ntd++;
04293 }
04294 else
04295 {
04296 pv[i].setSide(-2);
04297 rme0 += pv[i].getMass();
04298 nte++;
04299 }
04300 }
04301 else
04302 {
04303 rme0 += pv[i].getMass();
04304 nte++;
04305 }
04306 }
04307
04308 G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
04309 G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
04310
04311 G4double rmc = rmc0, rmd = rmd0, rme = rme0;
04312 G4int ntc1 = Imin(4,ntc-1);
04313 G4int ntd1 = Imin(4,ntd-1);
04314 G4int nte1 = Imin(4,nte-1);
04315 if (ntc > 1) rmc = rmc0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntc1])/gpar[ntc1];
04316 if (ntd > 1) rmd = rmd0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntd1])/gpar[ntd1];
04317 if (nte > 1) rme = rme0 + std::pow(-std::log(1.-G4UniformRand()),cpar[nte1])/gpar[nte1];
04318 while( (rmc+rmd) > centerOfMassEnergy)
04319 {
04320 if ((rmc == rmc0) && (rmd == rmd0))
04321 {
04322 rmd *= 0.999*centerOfMassEnergy/(rmc+rmd);
04323 rmc *= 0.999*centerOfMassEnergy/(rmc+rmd);
04324 }
04325 else
04326 {
04327 rmc = 0.1*rmc0 + 0.9*rmc;
04328 rmd = 0.1*rmd0 + 0.9*rmd;
04329 }
04330 }
04331 if(verboseLevel > 1)
04332 G4cout << " Cluster Masses: " << ntc << " " << rmc << " " << ntd << " "
04333 << rmd << " " << nte << " " << rme << G4endl;
04334
04335
04336 G4HEVector* pvmx = new G4HEVector[11];
04337
04338 pvmx[1].setMass( incidentMass);
04339 pvmx[1].setMomentumAndUpdate( 0., 0., incidentTotalMomentum);
04340 pvmx[2].setMass( targetMass);
04341 pvmx[2].setMomentumAndUpdate( 0., 0., 0.);
04342 pvmx[0].Add( pvmx[1], pvmx[2] );
04343 pvmx[1].Lor( pvmx[1], pvmx[0] );
04344 pvmx[2].Lor( pvmx[2], pvmx[0] );
04345
04346 G4double pf = std::sqrt(Amax(0.0001, sqr(sqr(centerOfMassEnergy) + rmd*rmd -rmc*rmc)
04347 - 4*sqr(centerOfMassEnergy)*rmd*rmd))/(2.*centerOfMassEnergy);
04348 pvmx[3].setMass( rmc );
04349 pvmx[4].setMass( rmd );
04350 pvmx[3].setEnergy( std::sqrt(pf*pf + rmc*rmc) );
04351 pvmx[4].setEnergy( std::sqrt(pf*pf + rmd*rmd) );
04352
04353 G4double tvalue = -DBL_MAX;
04354 G4double bvalue = Amax(0.01, 4.0 + 1.6*std::log(incidentTotalMomentum));
04355 if (bvalue != 0.0) tvalue = std::log(G4UniformRand())/bvalue;
04356 G4double pin = pvmx[1].Length();
04357 G4double tacmin = sqr( pvmx[1].getEnergy() - pvmx[3].getEnergy()) - sqr( pin - pf);
04358 G4double ctet = Amax(-1., Amin(1., 1.+2.*(tvalue-tacmin)/Amax(1.e-10, 4.*pin*pf)));
04359 G4double stet = std::sqrt(Amax(0., 1.0 - ctet*ctet));
04360 G4double phi = twopi * G4UniformRand();
04361 pvmx[3].setMomentum( pf * stet * std::sin(phi),
04362 pf * stet * std::cos(phi),
04363 pf * ctet );
04364 pvmx[4].Smul( pvmx[3], -1.);
04365
04366 if (nte > 0)
04367 {
04368 G4double ekit1 = 0.04;
04369 G4double ekit2 = 0.6;
04370 G4double gaval = 1.2;
04371 if (incidentKineticEnergy <= 5.)
04372 {
04373 ekit1 *= sqr(incidentKineticEnergy)/25.;
04374 ekit2 *= sqr(incidentKineticEnergy)/25.;
04375 }
04376 G4double avalue = (1.-gaval)/(std::pow(ekit2, 1.-gaval)-std::pow(ekit1, 1.-gaval));
04377 for (i=0; i < vecLen; i++)
04378 {
04379 if (pv[i].getSide() == -2)
04380 {
04381 G4double ekit = std::pow(G4UniformRand()*(1.-gaval)/avalue +std::pow(ekit1, 1.-gaval),
04382 1./(1.-gaval));
04383 pv[i].setKineticEnergyAndUpdate( ekit );
04384 ctet = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
04385 stet = std::sqrt( Amax( 0.0, 1. - ctet*ctet ));
04386 phi = G4UniformRand()*twopi;
04387 G4double pp = pv[i].Length();
04388 pv[i].setMomentum( pp * stet * std::sin(phi),
04389 pp * stet * std::cos(phi),
04390 pp * ctet );
04391 pv[i].Lor( pv[i], pvmx[0] );
04392 }
04393 }
04394 }
04395
04396
04397 pvmx[5].SmulAndUpdate( pvmx[3], -1.);
04398 pvmx[6].SmulAndUpdate( pvmx[4], -1.);
04399
04400 if (verboseLevel > 1) {
04401 G4cout << " General vectors before Phase space Generation " << G4endl;
04402 for (i=0; i<7; i++) pvmx[i].Print(i);
04403 }
04404
04405
04406 G4HEVector* tempV = new G4HEVector[18];
04407 G4bool constantCrossSection = true;
04408 G4double wgt;
04409 G4int npg;
04410
04411 if (ntc > 1)
04412 {
04413 npg = 0;
04414 for (i=0; i < vecLen; i++)
04415 {
04416 if (pv[i].getSide() > 0)
04417 {
04418 tempV[npg++] = pv[i];
04419 if(verboseLevel > 1) pv[i].Print(i);
04420 }
04421 }
04422 wgt = NBodyPhaseSpace( pvmx[3].getMass(), constantCrossSection, tempV, npg);
04423
04424 npg = 0;
04425 for (i=0; i < vecLen; i++)
04426 {
04427 if (pv[i].getSide() > 0)
04428 {
04429 pv[i].setMomentum( tempV[npg++].getMomentum());
04430 pv[i].SmulAndUpdate( pv[i], 1. );
04431 pv[i].Lor( pv[i], pvmx[5] );
04432 if(verboseLevel > 1) pv[i].Print(i);
04433 }
04434 }
04435 }
04436 else if(ntc == 1)
04437 {
04438 for(i=0; i<vecLen; i++)
04439 {
04440 if(pv[i].getSide() > 0) pv[i].setMomentumAndUpdate(pvmx[3].getMomentum());
04441 if(verboseLevel > 1) pv[i].Print(i);
04442 }
04443 }
04444 else
04445 {
04446 }
04447
04448 if (ntd > 1)
04449 {
04450 npg = 0;
04451 for (i=0; i < vecLen; i++)
04452 {
04453 if (pv[i].getSide() == -1)
04454 {
04455 tempV[npg++] = pv[i];
04456 if(verboseLevel > 1) pv[i].Print(i);
04457 }
04458 }
04459 wgt = NBodyPhaseSpace( pvmx[4].getMass(), constantCrossSection, tempV, npg);
04460
04461 npg = 0;
04462 for (i=0; i < vecLen; i++)
04463 {
04464 if (pv[i].getSide() == -1)
04465 {
04466 pv[i].setMomentum( tempV[npg++].getMomentum());
04467 pv[i].SmulAndUpdate( pv[i], 1.);
04468 pv[i].Lor( pv[i], pvmx[6] );
04469 if(verboseLevel > 1) pv[i].Print(i);
04470 }
04471 }
04472 }
04473 else if(ntd == 1)
04474 {
04475 for(i=0; i<vecLen; i++)
04476 {
04477 if(pv[i].getSide() == -1) pv[i].setMomentumAndUpdate(pvmx[4].getMomentum());
04478 if(verboseLevel > 1) pv[i].Print(i);
04479 }
04480 }
04481 else
04482 {
04483 }
04484
04485 if(verboseLevel > 1)
04486 {
04487 G4cout << " Vectors after PhaseSpace generation " << G4endl;
04488 for(i=0;i<vecLen; i++) pv[i].Print(i);
04489 }
04490
04491
04492
04493 targ = 0;
04494 for( i=0; i < vecLen; i++ )
04495 {
04496 if( pv[i].getType() == baryonType )targ++;
04497 if( pv[i].getType() == antiBaryonType )targ++;
04498 pv[i].Lor( pv[i], pvmx[2] );
04499 }
04500 if (targ <1) targ =1;
04501
04502 if(verboseLevel > 1) {
04503 G4cout << " Transformation in Lab- System " << G4endl;
04504 for(i=0; i<vecLen; i++) pv[i].Print(i);
04505 }
04506
04507
04508
04509
04510 G4double ekin, teta;
04511
04512 if (lead) {
04513 for (i = 0; i < vecLen; i++) {
04514 if (pv[i].getCode() == lead) {
04515
04516
04517
04518
04519 break;
04520 }
04521 }
04522
04523
04524
04525
04526
04527
04528
04529
04530
04531
04532
04533
04534
04535
04536
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546 }
04547
04548 pvmx[4].setMass( incidentMass);
04549 pvmx[4].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
04550
04551 G4double ekin0 = pvmx[4].getKineticEnergy();
04552
04553 pvmx[5].setMass ( protonMass * targ);
04554 pvmx[5].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
04555
04556 ekin = pvmx[4].getEnergy() + pvmx[5].getEnergy();
04557
04558 pvmx[6].Add( pvmx[4], pvmx[5] );
04559 pvmx[4].Lor( pvmx[4], pvmx[6] );
04560 pvmx[5].Lor( pvmx[5], pvmx[6] );
04561
04562 G4double tecm = pvmx[4].getEnergy() + pvmx[5].getEnergy();
04563
04564 pvmx[8].setZero();
04565
04566 G4double ekin1 = 0.0;
04567
04568 for( i=0; i < vecLen; i++ )
04569 {
04570 pvmx[8].Add( pvmx[8], pv[i] );
04571 ekin1 += pv[i].getKineticEnergy();
04572 ekin -= pv[i].getMass();
04573 }
04574
04575 if( vecLen > 1 && vecLen < 19 )
04576 {
04577 constantCrossSection = true;
04578 G4HEVector pw[18];
04579 for(i=0;i<vecLen;i++) pw[i] = pv[i];
04580 wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
04581 ekin = 0.0;
04582 for( i=0; i < vecLen; i++ )
04583 {
04584 pvmx[7].setMass( pw[i].getMass());
04585 pvmx[7].setMomentum( pw[i].getMomentum() );
04586 pvmx[7].SmulAndUpdate( pvmx[7], 1.);
04587 pvmx[7].Lor( pvmx[7], pvmx[5] );
04588 ekin += pvmx[7].getKineticEnergy();
04589 }
04590 teta = pvmx[8].Ang( pvmx[4] );
04591 if (verboseLevel > 1)
04592 G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
04593 << " " << ekin1 << " " << ekin << G4endl;
04594 }
04595
04596 if( ekin1 != 0.0 )
04597 {
04598 pvmx[7].setZero();
04599 wgt = ekin/ekin1;
04600 ekin1 = 0.;
04601 for( i=0; i < vecLen; i++ )
04602 {
04603 pvMass = pv[i].getMass();
04604 ekin = pv[i].getKineticEnergy() * wgt;
04605 pv[i].setKineticEnergyAndUpdate( ekin );
04606 ekin1 += ekin;
04607 pvmx[7].Add( pvmx[7], pv[i] );
04608 }
04609 teta = pvmx[7].Ang( pvmx[4] );
04610 if (verboseLevel > 1)
04611 G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
04612 << ekin1 << G4endl;
04613 }
04614
04615
04616
04617 G4double ry = G4UniformRand();
04618 G4double rz = G4UniformRand();
04619 G4double rx = twopi*rz;
04620 G4double a1 = std::sqrt(-2.0*std::log(ry));
04621 G4double rantarg1 = a1*std::cos(rx)*0.02*targ/G4double(vecLen);
04622 G4double rantarg2 = a1*std::sin(rx)*0.02*targ/G4double(vecLen);
04623
04624 for (i = 0; i < vecLen; i++)
04625 pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
04626 pv[i].getMomentum().y()+rantarg2 );
04627
04628 if (verboseLevel > 1) {
04629 pvmx[7].setZero();
04630 for (i = 0; i < vecLen; i++) pvmx[7].Add( pvmx[7], pv[i] );
04631 teta = pvmx[7].Ang( pvmx[4] );
04632 G4cout << " After smearing " << teta << G4endl;
04633 }
04634
04635
04636
04637
04638
04639
04640
04641
04642 G4double dekin = 0.0;
04643 G4int npions = 0;
04644 G4double ek1 = 0.0;
04645 G4double alekw, xxh;
04646 G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
04647 G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
04648 G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
04649
04650 for (i = 0; i < vecLen; i++) {
04651 pv[i].Defs1( pv[i], pvI );
04652 if (atomicWeight > 1.5) {
04653 ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa*( 1. + 0.5*normal()));
04654 alekw = std::log( incidentKineticEnergy );
04655 xxh = 1.;
04656 xxh = 1.;
04657 if (incidentCode == pionPlusCode || incidentCode == pionMinusCode) {
04658 if (pv[i].getCode() == pionZeroCode) {
04659 if (G4UniformRand() < std::log(atomicWeight)) {
04660 if (alekw > alem[0]) {
04661 G4int jmax = 1;
04662 for (j = 1; j < 8; j++) {
04663 if (alekw < alem[j]) {
04664 jmax = j;
04665 break;
04666 }
04667 }
04668 xxh = (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alekw
04669 + val0[jmax-1] - (val0[jmax]-val0[jmax-1])/(alem[jmax]-alem[jmax-1])*alem[jmax-1];
04670 xxh = 1. - xxh;
04671 }
04672 }
04673 }
04674 }
04675 dekin += ekin*(1.-xxh);
04676 ekin *= xxh;
04677 pv[i].setKineticEnergyAndUpdate( ekin );
04678 pvCode = pv[i].getCode();
04679 if ((pvCode == pionPlusCode) ||
04680 (pvCode == pionMinusCode) ||
04681 (pvCode == pionZeroCode)) {
04682 npions += 1;
04683 ek1 += ekin;
04684 }
04685 }
04686 }
04687
04688 if( (ek1 > 0.0) && (npions > 0) )
04689 {
04690 dekin = 1.+dekin/ek1;
04691 for (i = 0; i < vecLen; i++)
04692 {
04693 pvCode = pv[i].getCode();
04694 if((pvCode == pionPlusCode) || (pvCode == pionMinusCode) || (pvCode == pionZeroCode))
04695 {
04696 ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
04697 pv[i].setKineticEnergyAndUpdate( ekin );
04698 }
04699 }
04700 }
04701 if (verboseLevel > 1)
04702 { G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
04703 for (i=0; i<vecLen; i++) pv[i].Print(i);
04704 }
04705
04706
04707
04708
04709
04710 if (verboseLevel > 1)
04711 G4cout << " Evaporation " << atomicWeight << " "
04712 << excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
04713
04714 if( atomicWeight > 1.5 )
04715 {
04716
04717 G4double sprob, cost, sint, ekin2, ran, pp, eka;
04718 G4int spall(0), nbl(0);
04719
04720
04721 if( incidentKineticEnergy < 5.0 )
04722 sprob = 0.0;
04723 else
04724
04725 sprob = Amin(1., 0.000314*atomicWeight*std::log(incidentKineticEnergy-4.));
04726
04727
04728 if( excitationEnergyGNP >= 0.001 )
04729 {
04730
04731
04732
04733 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyGNP/
04734 (excitationEnergyGNP+excitationEnergyDTA));
04735 if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
04736 if (verboseLevel > 1)
04737 G4cout << " evaporation " << targ << " " << nbl << " "
04738 << sprob << G4endl;
04739 spall = targ;
04740 if( nbl > 0)
04741 {
04742 ekin = excitationEnergyGNP/nbl;
04743 ekin2 = 0.0;
04744 for( i=0; i<nbl; i++ )
04745 {
04746 if( G4UniformRand() < sprob ) continue;
04747 if( ekin2 > excitationEnergyGNP) break;
04748 ran = G4UniformRand();
04749 ekin1 = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
04750 if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
04751 ekin2 += ekin1;
04752 if( ekin2 > excitationEnergyGNP )
04753 ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
04754 if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
04755 pv[vecLen].setDefinition( "Proton");
04756 else
04757 pv[vecLen].setDefinition( "Neutron");
04758 spall++;
04759 cost = G4UniformRand() * 2.0 - 1.0;
04760 sint = std::sqrt(std::fabs(1.0-cost*cost));
04761 phi = twopi * G4UniformRand();
04762 pv[vecLen].setFlag( true );
04763 pv[vecLen].setSide( -4 );
04764 pvMass = pv[vecLen].getMass();
04765 pv[vecLen].setTOF( 1.0 );
04766 pvEnergy = ekin1 + pvMass;
04767 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
04768 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
04769 pp*sint*std::cos(phi),
04770 pp*cost );
04771 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
04772 vecLen++;
04773 }
04774 if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
04775 {
04776 G4int ika, kk = 0;
04777 eka = incidentKineticEnergy;
04778 if( eka > 1.0 )eka *= eka;
04779 eka = Amax( 0.1, eka );
04780 ika = G4int(3.6*std::exp((atomicNumber*atomicNumber
04781 /atomicWeight-35.56)/6.45)/eka);
04782 if( ika > 0 )
04783 {
04784 for( i=(vecLen-1); i>=0; i-- )
04785 {
04786 if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
04787 {
04788 G4HEVector pTemp = pv[i];
04789 pv[i].setDefinition( "Neutron");
04790 pv[i].setMomentumAndUpdate(pTemp.getMomentum());
04791 if (verboseLevel > 1) pv[i].Print(i);
04792 if( ++kk > ika ) break;
04793 }
04794 }
04795 }
04796 }
04797 }
04798 }
04799
04800
04801
04802
04803 if( excitationEnergyDTA >= 0.001 )
04804 {
04805 nbl = Poisson( (1.5+1.25*targ)*excitationEnergyDTA
04806 /(excitationEnergyGNP+excitationEnergyDTA));
04807
04808
04809
04810 if( nbl > 0 )
04811 {
04812 ekin = excitationEnergyDTA/nbl;
04813 ekin2 = 0.0;
04814 for( i=0; i<nbl; i++ )
04815 {
04816 if( G4UniformRand() < sprob ) continue;
04817 if( ekin2 > excitationEnergyDTA) break;
04818 ran = G4UniformRand();
04819 ekin1 = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
04820 if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
04821 ekin2 += ekin1;
04822 if( ekin2 > excitationEnergyDTA)
04823 ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
04824 cost = G4UniformRand()*2.0 - 1.0;
04825 sint = std::sqrt(std::fabs(1.0-cost*cost));
04826 phi = twopi*G4UniformRand();
04827 ran = G4UniformRand();
04828 if( ran <= 0.60 )
04829 pv[vecLen].setDefinition( "Deuteron");
04830 else if (ran <= 0.90)
04831 pv[vecLen].setDefinition( "Triton");
04832 else
04833 pv[vecLen].setDefinition( "Alpha");
04834 spall += (int)(pv[vecLen].getMass() * 1.066);
04835 if( spall > atomicWeight ) break;
04836 pv[vecLen].setFlag( true );
04837 pv[vecLen].setSide( -4 );
04838 pvMass = pv[vecLen].getMass();
04839 pv[vecLen].setTOF( 1.0 );
04840 pvEnergy = pvMass + ekin1;
04841 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
04842 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
04843 pp*sint*std::cos(phi),
04844 pp*cost );
04845 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
04846 vecLen++;
04847 }
04848 }
04849 }
04850 }
04851 if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
04852 {
04853 for( i=0; i<vecLen; i++ )
04854 {
04855 G4double etb = pv[i].getKineticEnergy();
04856 if( etb >= incidentKineticEnergy )
04857 pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
04858 }
04859 }
04860
04861
04862
04863 G4double tof = incidentTOF;
04864 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
04865 && (incidentKineticEnergy <= 0.2) )
04866 tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
04867 for ( i=0; i < vecLen; i++)
04868 {
04869
04870 pv[i].setTOF ( tof );
04871
04872 }
04873
04874 for(i=0; i<vecLen; i++)
04875 {
04876 if(pv[i].getName() == "KaonZero" || pv[i].getName() == "AntiKaonZero")
04877 {
04878 pvmx[0] = pv[i];
04879 if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
04880 else pv[i].setDefinition("KaonZeroLong");
04881 pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
04882 }
04883 }
04884
04885 successful = true;
04886 delete [] pvmx;
04887 delete [] tempV;
04888 return;
04889 }
04890
04891 void
04892 G4HEInelastic::QuasiElasticScattering(G4bool& successful,
04893 G4HEVector pv[],
04894 G4int& vecLen,
04895 G4double& excitationEnergyGNP,
04896 G4double& excitationEnergyDTA,
04897 const G4HEVector& incidentParticle,
04898 const G4HEVector& targetParticle,
04899 G4double atomicWeight,
04900 G4double atomicNumber)
04901 {
04902
04903
04904
04905
04906
04907 G4int protonCode = Proton.getCode();
04908 G4String mesonType = PionPlus.getType();
04909 G4String baryonType = Proton.getType();
04910 G4String antiBaryonType = AntiProton.getType();
04911
04912 G4double targetMass = targetParticle.getMass();
04913 G4double incidentMass = incidentParticle.getMass();
04914 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
04915 G4double incidentEnergy = incidentParticle.getEnergy();
04916 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
04917 G4String incidentType = incidentParticle.getType();
04918
04919 G4double incidentTOF = 0.;
04920
04921
04922 G4int i;
04923
04924 if (verboseLevel > 1)
04925 G4cout << " G4HEInelastic::QuasiElasticScattering " << G4endl;
04926
04927 if (incidentTotalMomentum < 0.01 || vecLen < 2) {
04928 successful = false;
04929 return;
04930 }
04931
04932 G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
04933 +2.*targetMass*incidentEnergy);
04934
04935 G4HEVector pvI = incidentParticle;
04936 pvI.setSide( 1 );
04937
04938 G4HEVector pvT = targetParticle;
04939 pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
04940 pvT.setSide( -1 );
04941 pvT.setTOF( -1.);
04942
04943 G4HEVector* pvmx = new G4HEVector[3];
04944
04945 if (atomicWeight > 1.5) {
04946
04947 if ( (pvI.getCode() == pv[0].getCode() )
04948 && (pvT.getCode() == pv[1].getCode() )
04949 && (excitationEnergyGNP < 0.001)
04950 && (excitationEnergyDTA < 0.001) ) {
04951 successful = false;
04952 delete [] pvmx;
04953 return;
04954 }
04955 }
04956
04957 pv[0].setSide( 1 );
04958 pv[0].setFlag( false );
04959 pv[0].setTOF( incidentTOF);
04960 pv[0].setMomentumAndUpdate( incidentParticle.getMomentum() );
04961 pv[1].setSide( -1 );
04962 pv[1].setFlag( false );
04963 pv[1].setTOF( incidentTOF);
04964 pv[1].setMomentumAndUpdate(targetParticle.getMomentum() );
04965
04966 if ((incidentTotalMomentum > 0.1) && (centerOfMassEnergy > 0.01) ) {
04967 if (pv[1].getType() == mesonType) {
04968 if (G4UniformRand() < 0.5)
04969 pv[1].setDefinition( "Proton");
04970 else
04971 pv[1].setDefinition( "Neutron");
04972 }
04973 pvmx[0].Add( pvI, pvT );
04974 pvmx[1].Lor( pvI, pvmx[0] );
04975 pvmx[2].Lor( pvT, pvmx[0] );
04976 G4double pin = pvmx[1].Length();
04977 G4double bvalue = Amax(0.01 , 4.225+1.795*std::log(incidentTotalMomentum));
04978 G4double pf = sqr(sqr(centerOfMassEnergy) + sqr(pv[1].getMass()) - sqr(pv[0].getMass()))
04979 - 4 * sqr(centerOfMassEnergy) * sqr(pv[1].getMass());
04980
04981 if (pf < 0.001) {
04982 successful = false;
04983 delete [] pvmx;
04984 return;
04985 }
04986 pf = std::sqrt(pf)/(2.*centerOfMassEnergy);
04987 G4double btrang = 4. * bvalue * pin * pf;
04988 G4double exindt = -1.;
04989 if (btrang < 46.) exindt += std::exp(-btrang);
04990 G4double tdn = std::log(1. + G4UniformRand()*exindt)/btrang;
04991 G4double ctet = Amax( -1., Amin(1., 1. + 2.*tdn));
04992 G4double stet = std::sqrt((1.-ctet)*(1.+ctet));
04993 G4double phi = twopi * G4UniformRand();
04994 pv[0].setMomentumAndUpdate( pf*stet*std::sin(phi),
04995 pf*stet*std::cos(phi),
04996 pf*ctet );
04997 pv[1].SmulAndUpdate( pv[0], -1.);
04998 for (i = 0; i < 2; i++) {
04999
05000
05001 pv[i].Lor(pv[i], pvmx[0]);
05002 pv[i].Defs1(pv[i], pvI);
05003 if (atomicWeight > 1.5) {
05004 G4double ekin = pv[i].getKineticEnergy()
05005 - 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.)
05006 *(1. + 0.5*normal());
05007 ekin = Amax(0.0001, ekin);
05008 pv[i].setKineticEnergyAndUpdate( ekin );
05009 }
05010 }
05011 }
05012 vecLen = 2;
05013
05014
05015
05016
05017
05018 if (verboseLevel > 1)
05019 G4cout << " Evaporation " << atomicWeight << " "
05020 << excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
05021
05022 if( atomicWeight > 1.5 )
05023 {
05024
05025 G4double sprob, cost, sint, ekin2, ran, pp, eka;
05026 G4double ekin, cfa, ekin1, phi, pvMass, pvEnergy;
05027 G4int spall(0), nbl(0);
05028
05029
05030 sprob = 0.;
05031 cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
05032
05033
05034 if( excitationEnergyGNP >= 0.001 )
05035 {
05036
05037
05038
05039 nbl = Poisson( excitationEnergyGNP/0.02);
05040 if( nbl > atomicWeight ) nbl = (int)(atomicWeight);
05041 if (verboseLevel > 1)
05042 G4cout << " evaporation " << nbl << " " << sprob << G4endl;
05043 spall = 0;
05044 if( nbl > 0)
05045 {
05046 ekin = excitationEnergyGNP/nbl;
05047 ekin2 = 0.0;
05048 for( i=0; i<nbl; i++ )
05049 {
05050 if( G4UniformRand() < sprob ) continue;
05051 if( ekin2 > excitationEnergyGNP) break;
05052 ran = G4UniformRand();
05053 ekin1 = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
05054 if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
05055 ekin2 += ekin1;
05056 if( ekin2 > excitationEnergyGNP)
05057 ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
05058 if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
05059 pv[vecLen].setDefinition( "Proton");
05060 else
05061 pv[vecLen].setDefinition( "Neutron");
05062 spall++;
05063 cost = G4UniformRand() * 2.0 - 1.0;
05064 sint = std::sqrt(std::fabs(1.0-cost*cost));
05065 phi = twopi * G4UniformRand();
05066 pv[vecLen].setFlag( true );
05067 pv[vecLen].setSide( -4 );
05068 pvMass = pv[vecLen].getMass();
05069 pv[vecLen].setTOF( 1.0 );
05070 pvEnergy = ekin1 + pvMass;
05071 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
05072 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
05073 pp*sint*std::cos(phi),
05074 pp*cost );
05075 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
05076 vecLen++;
05077 }
05078 if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
05079 {
05080 G4int ika, kk = 0;
05081 eka = incidentKineticEnergy;
05082 if( eka > 1.0 )eka *= eka;
05083 eka = Amax( 0.1, eka );
05084 ika = G4int(3.6*std::exp((atomicNumber*atomicNumber
05085 /atomicWeight-35.56)/6.45)/eka);
05086 if( ika > 0 )
05087 {
05088 for( i=(vecLen-1); i>=0; i-- )
05089 {
05090 if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
05091 {
05092 pv[i].setDefinition( "Neutron" );
05093 if (verboseLevel > 1) pv[i].Print(i);
05094 if( ++kk > ika ) break;
05095 }
05096 }
05097 }
05098 }
05099 }
05100 }
05101
05102
05103
05104
05105 if( excitationEnergyDTA >= 0.001 )
05106 {
05107 nbl = (G4int)(2.*std::log(atomicWeight));
05108
05109
05110
05111 if( nbl > 0 )
05112 {
05113 ekin = excitationEnergyDTA/nbl;
05114 ekin2 = 0.0;
05115 for( i=0; i<nbl; i++ )
05116 {
05117 if( G4UniformRand() < sprob ) continue;
05118 if( ekin2 > excitationEnergyDTA) break;
05119 ran = G4UniformRand();
05120 ekin1 = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
05121 if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
05122 ekin2 += ekin1;
05123 if( ekin2 > excitationEnergyDTA)
05124 ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
05125 cost = G4UniformRand()*2.0 - 1.0;
05126 sint = std::sqrt(std::fabs(1.0-cost*cost));
05127 phi = twopi*G4UniformRand();
05128 ran = G4UniformRand();
05129 if( ran <= 0.60 )
05130 pv[vecLen].setDefinition( "Deuteron");
05131 else if (ran <= 0.90)
05132 pv[vecLen].setDefinition( "Triton");
05133 else
05134 pv[vecLen].setDefinition( "Alpha");
05135 spall += (int)(pv[vecLen].getMass() * 1.066);
05136 if( spall > atomicWeight ) break;
05137 pv[vecLen].setFlag( true );
05138 pv[vecLen].setSide( -4 );
05139 pvMass = pv[vecLen].getMass();
05140 pv[vecLen].setTOF( 1.0 );
05141 pvEnergy = pvMass + ekin1;
05142 pp = std::sqrt( std::fabs( pvEnergy*pvEnergy - pvMass*pvMass ) );
05143 pv[vecLen].setMomentumAndUpdate( pp*sint*std::sin(phi),
05144 pp*sint*std::cos(phi),
05145 pp*cost );
05146 if (verboseLevel > 1) pv[vecLen].Print(vecLen);
05147 vecLen++;
05148 }
05149 }
05150 }
05151 }
05152
05153
05154
05155 G4double tof = incidentTOF;
05156 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
05157 && (incidentKineticEnergy <= 0.2) )
05158 tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
05159 for ( i=0; i < vecLen; i++)
05160 {
05161
05162 pv[i].setTOF ( tof );
05163
05164 }
05165
05166 for(i=0; i<vecLen; i++)
05167 {
05168 if(pv[i].getName() == "KaonZero" || pv[i].getName() == "AntiKaonZero")
05169 {
05170 pvmx[0] = pv[i];
05171 if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
05172 else pv[i].setDefinition("KaonZeroLong");
05173 pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
05174 }
05175 }
05176
05177 successful = true;
05178 delete [] pvmx;
05179 return;
05180 }
05181
05182 void
05183 G4HEInelastic::ElasticScattering(G4bool& successful,
05184 G4HEVector pv[],
05185 G4int& vecLen,
05186 const G4HEVector& incidentParticle,
05187 G4double atomicWeight,
05188 G4double )
05189 {
05190 if (verboseLevel > 1)
05191 G4cout << " G4HEInelastic::ElasticScattering " << G4endl;
05192
05193 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
05194 if (verboseLevel > 1)
05195 G4cout << "DoIt: Incident particle momentum="
05196 << incidentTotalMomentum << " GeV" << G4endl;
05197 if (incidentTotalMomentum < 0.01) {
05198 successful = false;
05199 return;
05200 }
05201
05202 if (atomicWeight < 0.5)
05203 {
05204 successful = false;
05205 return;
05206 }
05207 pv[0] = incidentParticle;
05208 vecLen = 1;
05209
05210 G4double aa, bb, cc, dd, rr;
05211 if (atomicWeight <= 62.)
05212 {
05213 aa = std::pow(atomicWeight, 1.63);
05214 bb = 14.5*std::pow(atomicWeight, 0.66);
05215 cc = 1.4*std::pow(atomicWeight, 0.33);
05216 dd = 10.;
05217 }
05218 else
05219 {
05220 aa = std::pow(atomicWeight, 1.33);
05221 bb = 60.*std::pow(atomicWeight, 0.33);
05222 cc = 0.4*std::pow(atomicWeight, 0.40);
05223 dd = 10.;
05224 }
05225 aa = aa/bb;
05226 cc = cc/dd;
05227 G4double ran = G4UniformRand();
05228 rr = (aa + cc)*ran;
05229 if (verboseLevel > 1)
05230 {
05231 G4cout << "ElasticScattering: aa,bb,cc,dd,rr" << G4endl;
05232 G4cout << aa << " " << bb << " " << cc << " " << dd << " "
05233 << rr << G4endl;
05234 }
05235 G4double t1 = -std::log(ran)/bb;
05236 G4double t2 = -std::log(ran)/dd;
05237 if (verboseLevel > 1) {
05238 G4cout << "t1,fctcos " << t1 << " " << fctcos(t1, aa, bb, cc, dd, rr)
05239 << G4endl;
05240 G4cout << "t2,fctcos " << t2 << " " << fctcos(t2, aa, bb, cc, dd, rr)
05241 << G4endl;
05242 }
05243 G4double eps = 0.001;
05244 G4int ind1 = 10;
05245 G4double t;
05246 G4int ier1;
05247 ier1 = rtmi(&t, t1, t2, eps, ind1, aa, bb, cc, dd, rr);
05248 if (verboseLevel > 1) {
05249 G4cout << "From rtmi, ier1=" << ier1 << G4endl;
05250 G4cout << "t, fctcos " << t << " " << fctcos(t, aa, bb, cc, dd, rr)
05251 << G4endl;
05252 }
05253 if (ier1 != 0) t = 0.25*(3.*t1 + t2);
05254 if (verboseLevel > 1)
05255 G4cout << "t, fctcos " << t << " " << fctcos(t, aa, bb, cc, dd, rr)
05256 << G4endl;
05257
05258 G4double phi = G4UniformRand()*twopi;
05259 rr = 0.5*t/sqr(incidentTotalMomentum);
05260 if (rr > 1.) rr = 0.;
05261 if (verboseLevel > 1)
05262 G4cout << "rr=" << rr << G4endl;
05263 G4double cost = 1. - rr;
05264 G4double sint = std::sqrt(Amax(rr*(2. - rr), 0.));
05265 if (verboseLevel > 1)
05266 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
05267
05268 G4HEVector pv0;
05269 G4HEVector pvI;
05270 pvI.setMass( incidentParticle.getMass() );
05271 pvI.setMomentum( incidentParticle.getMomentum() );
05272 pvI.SmulAndUpdate( pvI, 1. );
05273 pv0.setMass( pvI.getMass() );
05274
05275 pv0.setMomentumAndUpdate( incidentTotalMomentum * sint * std::sin(phi),
05276 incidentTotalMomentum * sint * std::cos(phi),
05277 incidentTotalMomentum * cost );
05278 pv0.Defs1( pv0, pvI );
05279
05280 successful = true;
05281 return;
05282 }
05283
05284
05285 G4int
05286 G4HEInelastic::rtmi(G4double *x, G4double xli, G4double xri, G4double eps,
05287 G4int iend,
05288 G4double aa, G4double bb, G4double cc, G4double dd,
05289 G4double rr)
05290 {
05291 G4int ier = 0;
05292 G4double xl = xli;
05293 G4double xr = xri;
05294 *x = xl;
05295 G4double tol = *x;
05296 G4double f = fctcos(tol, aa, bb, cc, dd, rr);
05297 if (f == 0.) return ier;
05298 G4double fl, fr;
05299 fl = f;
05300 *x = xr;
05301 tol = *x;
05302 f = fctcos(tol, aa, bb, cc, dd, rr);
05303 if (f == 0.) return ier;
05304 fr = f;
05305
05306
05307 if (fl*fr >= 0.)
05308 {
05309 ier = 2;
05310 return ier;
05311 }
05312
05313
05314
05315 G4int i = 0;
05316 G4double tolf = 100.*eps;
05317
05318
05319
05320 label4:
05321 i++;
05322
05323
05324
05325 for (G4int k = 1; k <= iend; k++)
05326 {
05327 *x = 0.5*(xl + xr);
05328 tol = *x;
05329 f = fctcos(tol, aa, bb, cc, dd, rr);
05330 if (f == 0.) return 0;
05331 if (f*fr < 0.)
05332 {
05333 tol = xl;
05334 xl = xr;
05335 xr = tol;
05336 tol = fl;
05337 fl = fr;
05338 fr = tol;
05339 }
05340 tol = f - fl;
05341 G4double a = f*tol;
05342 a = a + a;
05343 if (a < fr*(fr - fl) && i <= iend) goto label17;
05344 xr = *x;
05345 fr = f;
05346
05347
05348 tol = eps;
05349 a = std::fabs(xr);
05350 if (a > 1.) tol = tol*a;
05351 if (std::fabs(xr - xl) <= tol && std::fabs(fr - fl) <= tolf) goto label14;
05352 }
05353
05354
05355
05356
05357
05358 ier = 1;
05359
05360 label14:
05361 if (std::fabs(fr) > std::fabs(fl))
05362 {
05363 *x = xl;
05364 f = fl;
05365 }
05366 return ier;
05367
05368
05369 label17:
05370 G4double a = fr - f;
05371 G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
05372 G4double xm = *x;
05373 G4double fm = f;
05374 *x = xl - dx;
05375 tol = *x;
05376 f = fctcos(tol, aa, bb, cc, dd, rr);
05377 if (f == 0.) return ier;
05378
05379
05380 tol = eps;
05381 a = std::fabs(*x);
05382 if (a > 1) tol = tol*a;
05383 if (std::fabs(dx) <= tol && std::fabs(f) <= tolf) return ier;
05384
05385
05386 if (f*fl < 0.)
05387 {
05388 xr = *x;
05389 fr = f;
05390 }
05391 else
05392 {
05393 xl = *x;
05394 fl = f;
05395 xr = xm;
05396 fr = fm;
05397 }
05398 goto label4;
05399 }
05400
05401
05402
05403
05404 G4double
05405 G4HEInelastic::fctcos(G4double t, G4double aa, G4double bb, G4double cc,
05406 G4double dd, G4double rr)
05407 {
05408 const G4double expxl = -82.;
05409 const G4double expxu = 82.;
05410
05411 G4double test1 = -bb*t;
05412 if (test1 > expxu) test1 = expxu;
05413 if (test1 < expxl) test1 = expxl;
05414
05415 G4double test2 = -dd*t;
05416 if (test2 > expxu) test2 = expxu;
05417 if (test2 < expxl) test2 = expxl;
05418
05419 return aa*std::exp(test1) + cc*std::exp(test2) - rr;
05420 }
05421
05422 G4double
05423 G4HEInelastic::NBodyPhaseSpace(const G4double totalEnergy,
05424 const G4bool constantCrossSection,
05425 G4HEVector vec[],
05426 G4int& vecLen )
05427 {
05428
05429
05430 G4int i;
05431
05432 const G4double expxu = std::log(FLT_MAX);
05433 const G4double expxl = -expxu;
05434
05435 if( vecLen < 2 ) {
05436 G4cerr << "*** Error in G4HEInelastic::GenerateNBodyEvent" << G4endl;
05437 G4cerr << " number of particles < 2" << G4endl;
05438 G4cerr << "totalEnergy = " << totalEnergy << ", vecLen = "
05439 << vecLen << G4endl;
05440 return -1.0;
05441 }
05442
05443 G4double* mass = new G4double [vecLen];
05444 G4double* energy = new G4double [vecLen];
05445 G4double** pcm;
05446 pcm = new G4double* [3];
05447 for( i=0; i<3; ++i )pcm[i] = new G4double [vecLen];
05448
05449 G4double totalMass = 0.0;
05450 G4double* sm = new G4double [vecLen];
05451
05452 for( i=0; i<vecLen; ++i ) {
05453 mass[i] = vec[i].getMass();
05454 vec[i].setMomentum( 0.0, 0.0, 0.0 );
05455 pcm[0][i] = 0.0;
05456 pcm[1][i] = 0.0;
05457 pcm[2][i] = 0.0;
05458 energy[i] = mass[i];
05459 totalMass += mass[i];
05460 sm[i] = totalMass;
05461 }
05462
05463 if (totalMass >= totalEnergy ) {
05464 if (verboseLevel > 1) {
05465 G4cout << "*** Error in G4HEInelastic::GenerateNBodyEvent" << G4endl;
05466 G4cout << " total mass (" << totalMass << ") >= total energy ("
05467 << totalEnergy << ")" << G4endl;
05468 }
05469 delete [] mass;
05470 delete [] energy;
05471 for( i=0; i<3; ++i )delete [] pcm[i];
05472 delete [] pcm;
05473 delete [] sm;
05474 return -1.0;
05475 }
05476
05477 G4double kineticEnergy = totalEnergy - totalMass;
05478 G4double* emm = new G4double [vecLen];
05479 emm[0] = mass[0];
05480 if (vecLen > 3) {
05481 G4double* ran = new G4double [vecLen];
05482 for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
05483 for( i=0; i<vecLen-1; ++i ) {
05484 for( G4int j=vecLen-1; j > i; --j ) {
05485 if( ran[i] > ran[j] ) {
05486 G4double temp = ran[i];
05487 ran[i] = ran[j];
05488 ran[j] = temp;
05489 }
05490 }
05491 }
05492 for( i=1; i<vecLen; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
05493 delete [] ran;
05494 } else {
05495 emm[1] = G4UniformRand()*kineticEnergy + sm[1];
05496 }
05497
05498 emm[vecLen-1] = totalEnergy;
05499
05500
05501
05502 G4bool lzero = true;
05503 G4double wtmax = 0.0;
05504 if (constantCrossSection) {
05505 G4double emmax = kineticEnergy + mass[0];
05506 G4double emmin = 0.0;
05507 for( i=1; i<vecLen; ++i ) {
05508 emmin += mass[i-1];
05509 emmax += mass[i];
05510 G4double wtfc = 0.0;
05511 if( emmax*emmax > 0.0 ) {
05512 G4double arg = emmax*emmax
05513 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
05514 - 2.0*(emmin*emmin+mass[i]*mass[i]);
05515 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
05516 }
05517 if( wtfc == 0.0 ) {
05518 lzero = false;
05519 break;
05520 }
05521 wtmax += std::log( wtfc );
05522 }
05523 if( lzero )
05524 wtmax = -wtmax;
05525 else
05526 wtmax = expxu;
05527 } else {
05528 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) *
05529 pi * std::pow( twopi, vecLen-2 ) / Factorial(vecLen-2) );
05530 }
05531 lzero = true;
05532 G4double* pd = new G4double [vecLen-1];
05533 for( i=0; i<vecLen-1; ++i ) {
05534 pd[i] = 0.0;
05535 if( emm[i+1]*emm[i+1] > 0.0 ) {
05536 G4double arg = emm[i+1]*emm[i+1]
05537 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
05538 /(emm[i+1]*emm[i+1])
05539 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
05540 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
05541 }
05542 if( pd[i] == 0.0 )
05543 lzero = false;
05544 else
05545 wtmax += std::log( pd[i] );
05546 }
05547 G4double weight = 0.0;
05548 if( lzero )weight = std::exp( Amax(Amin(wtmax,expxu),expxl) );
05549
05550 G4double bang, cb, sb, s0, s1, s2, c, esys, a, b, gama, beta;
05551 G4double ss;
05552 pcm[0][0] = 0.0;
05553 pcm[1][0] = pd[0];
05554 pcm[2][0] = 0.0;
05555 for( i=1; i<vecLen; ++i ) {
05556 pcm[0][i] = 0.0;
05557 pcm[1][i] = -pd[i-1];
05558 pcm[2][i] = 0.0;
05559 bang = twopi*G4UniformRand();
05560 cb = std::cos(bang);
05561 sb = std::sin(bang);
05562 c = 2.0*G4UniformRand() - 1.0;
05563 ss = std::sqrt( std::fabs( 1.0-c*c ) );
05564 if( i < vecLen-1 ) {
05565 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
05566 beta = pd[i]/esys;
05567 gama = esys/emm[i];
05568 for( G4int j=0; j<=i; ++j ) {
05569 s0 = pcm[0][j];
05570 s1 = pcm[1][j];
05571 s2 = pcm[2][j];
05572 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
05573 a = s0*c - s1*ss;
05574 pcm[1][j] = s0*ss + s1*c;
05575 b = pcm[2][j];
05576 pcm[0][j] = a*cb - b*sb;
05577 pcm[2][j] = a*sb + b*cb;
05578 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
05579 }
05580 } else {
05581 for( G4int j=0; j<=i; ++j ) {
05582 s0 = pcm[0][j];
05583 s1 = pcm[1][j];
05584 s2 = pcm[2][j];
05585 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
05586 a = s0*c - s1*ss;
05587 pcm[1][j] = s0*ss + s1*c;
05588 b = pcm[2][j];
05589 pcm[0][j] = a*cb - b*sb;
05590 pcm[2][j] = a*sb + b*cb;
05591 }
05592 }
05593 }
05594
05595 G4double pModule;
05596 for (i = 0; i < vecLen; ++i) {
05597 kineticEnergy = energy[i] - mass[i];
05598 pModule = std::sqrt( sqr(kineticEnergy) + 2*kineticEnergy*mass[i] );
05599 vec[i].setMomentum(pcm[0][i]/pModule,
05600 pcm[1][i]/pModule,
05601 pcm[2][i]/pModule);
05602 vec[i].setKineticEnergyAndUpdate( kineticEnergy );
05603 }
05604
05605 delete [] mass;
05606 delete [] energy;
05607 for( i=0; i<3; ++i )delete [] pcm[i];
05608 delete [] pcm;
05609 delete [] emm;
05610 delete [] sm;
05611 delete [] pd;
05612 return weight;
05613 }
05614
05615
05616 G4double
05617 G4HEInelastic::gpdk( G4double a, G4double b, G4double c )
05618 {
05619 if( a == 0.0 )
05620 {
05621 return 0.0;
05622 }
05623 else
05624 {
05625 G4double arg = a*a+(b*b-c*c)*(b*b-c*c)/(a*a)-2.0*(b*b+c*c);
05626 if( arg <= 0.0 )
05627 {
05628 return 0.0;
05629 }
05630 else
05631 {
05632 return 0.5*std::sqrt(std::fabs(arg));
05633 }
05634 }
05635 }
05636
05637
05638 G4double
05639 G4HEInelastic::NBodyPhaseSpace(G4int npart, G4HEVector pv[],
05640 G4double wmax, G4double wfcn,
05641 G4int maxtrial, G4int ntrial)
05642 { ntrial = 0;
05643 G4double wps(0);
05644 while ( ntrial < maxtrial)
05645 { ntrial++;
05646 G4int i, j;
05647 G4int nrn = 3*(npart-2)-4;
05648 G4double *ranarr = new G4double[nrn];
05649 for (i=0;i<nrn;i++) ranarr[i]=G4UniformRand();
05650 G4int nrnp = npart-4;
05651 if(nrnp > 1) QuickSort( ranarr, 0 , nrnp-1 );
05652 G4HEVector pvcms;
05653 pvcms.Add(pv[0],pv[1]);
05654 pvcms.Smul( pvcms, -1.);
05655 G4double rm = 0.;
05656 for (i=2;i<npart;i++) rm += pv[i].getMass();
05657 G4double rm1 = pvcms.getMass() - rm;
05658 rm -= pv[2].getMass();
05659 wps = (npart-3)*std::pow(rm1/sqr(twopi), npart-4)/(4*pi*pvcms.getMass());
05660 for (i=3; (i=npart-1);i++) wps /= i-2;
05661 G4double xxx = rm1/sqr(twopi);
05662 for (i=1; (i=npart-4); i++) wps /= xxx/i;
05663 wps /= (4*pi*pvcms.getMass());
05664 G4double p2,cost,sint,phi;
05665 j = 1;
05666 while (j)
05667 { j++;
05668 rm -= pv[j+1].getMass();
05669 if(j == npart-2) break;
05670 G4double rmass = rm + rm1*ranarr[npart-j-1];
05671 p2 = Alam(sqr(pvcms.getMass()), sqr(pv[j].getMass()),
05672 sqr(rmass))/(4.*sqr(pvcms.getMass()));
05673 cost = 1. - 2.*ranarr[npart+2*j-9];
05674 sint = std::sqrt(1.-cost*cost);
05675 phi = twopi*ranarr[npart+2*j-8];
05676 p2 = std::sqrt( Amax(0., p2));
05677 wps *= p2;
05678 pv[j].setMomentumAndUpdate( p2*sint*std::sin(phi), p2*sint*std::cos(phi),p2*cost);
05679 pv[j].Lor(pv[j], pvcms);
05680 pvcms.Add3( pvcms, pv[j] );
05681 pvcms.setEnergy(pvcms.getEnergy()-pv[j].getEnergy());
05682 pvcms.setMass( std::sqrt(sqr(pvcms.getEnergy()) - sqr(pvcms.Length())));
05683 }
05684 p2 = Alam(sqr(pvcms.getMass()), sqr(pv[j].getMass()),
05685 sqr(rm))/(4.*sqr(pvcms.getMass()));
05686 cost = 1. - 2.*ranarr[npart+2*j-9];
05687 sint = std::sqrt(1.-cost*cost);
05688 phi = twopi*ranarr[npart+2*j-8];
05689 p2 = std::sqrt( Amax(0. , p2));
05690 wps *= p2;
05691 pv[j].setMomentumAndUpdate( p2*sint*std::sin(phi), p2*sint*std::cos(phi), p2*cost);
05692 pv[j+1].setMomentumAndUpdate( -p2*sint*std::sin(phi), -p2*sint*std::cos(phi), -p2*cost);
05693 pv[j].Lor( pv[j], pvcms );
05694 pv[j+1].Lor( pv[j+1], pvcms );
05695 wfcn = CalculatePhaseSpaceWeight( npart );
05696 G4double wt = wps * wfcn;
05697 if (wt > wmax)
05698 { wmax = wt;
05699 G4cout << "maximum weight changed to " << wmax << G4endl;
05700 }
05701 wt = wt/wmax;
05702 if (G4UniformRand() < wt) break;
05703 }
05704 return wps;
05705 }
05706
05707
05708 void
05709 G4HEInelastic::QuickSort(G4double arr[], const G4int lidx, const G4int ridx)
05710 {
05711 G4double buffer;
05712 G4int k, e, mid;
05713 if(lidx>=ridx) return;
05714 mid = (int)((lidx+ridx)/2.);
05715 buffer = arr[lidx];
05716 arr[lidx]= arr[mid];
05717 arr[mid] = buffer;
05718 e = lidx;
05719 for (k=lidx+1;k<=ridx;k++)
05720 if (arr[k] < arr[lidx])
05721 { e++;
05722 buffer = arr[e];
05723 arr[e] = arr[k];
05724 arr[k] = buffer;
05725 }
05726 buffer = arr[lidx];
05727 arr[lidx]= arr[e];
05728 arr[e] = buffer;
05729 QuickSort(arr, lidx, e-1);
05730 QuickSort(arr, e+1 , ridx);
05731 return;
05732 }
05733
05734 G4double
05735 G4HEInelastic::Alam( G4double a, G4double b, G4double c)
05736 { return a*a + b*b + c*c - 2.*a*b - 2.*a*c -2.*b*c;
05737 }
05738
05739 G4double
05740 G4HEInelastic::CalculatePhaseSpaceWeight( G4int )
05741 { G4double wfcn = 1.;
05742 return wfcn;
05743 }
05744
05745 const std::pair<G4double, G4double> G4HEInelastic::GetFatalEnergyCheckLevels() const
05746 {
05747
05748 return std::pair<G4double, G4double>(5*perCent,250*GeV);
05749 }
05750
05751
05752
05753
05754
05755