124 for ( std::vector<G4PhysicsTable*>::iterator
it =
fAngleBank.begin();
127 if ( (*
it) ) (*it)->clearAndDestroy();
148 for( jEl = 0; jEl < numOfEl; ++jEl)
156 G4cout<<
"G4DiffuseElastic::Initialise() the element: "
157 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
206 if (iZ == 1 && iA == 1) theDef =
theProton;
209 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
210 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
224 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
226 if( cost >= 1.0 ) cost = 1.0;
227 else if( cost <= -1.0) cost = -1.0;
229 G4double thetaCMS = std::acos(cost);
261 if( z && (kRt > kRtC) )
293 if (iZ == 1 && iA == 1) theDef =
theProton;
296 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
297 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
311 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
313 if( cost >= 1.0 ) cost = 1.0;
314 else if( cost <= -1.0) cost = -1.0;
316 G4double thetaCMS = std::acos(cost);
343 if (iZ == 1 && iA == 1) theDef =
theProton;
346 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
347 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
361 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
363 if( cost >= 1.0 ) cost = 1.0;
364 else if( cost <= -1.0) cost = -1.0;
366 G4double thetaCMS = std::acos(cost);
386 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
396 diffuse = 0.63*
fermi;
404 diffuse = 0.63*
fermi;
415 diffuse = 0.63*
fermi;
426 bzero2 = bzero*bzero;
430 bonebyarg2 = bonebyarg*bonebyarg;
455 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
456 sigma += kr2*bonebyarg2;
474 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
487 bzero2 = bzero*bzero;
491 bonebyarg2 = bonebyarg*bonebyarg;
495 diffuse = 0.63*
fermi;
504 diffuse = 0.63*
fermi;
515 diffuse = 0.63*
fermi;
529 G4double sinHalfTheta = std::sin(0.5*theta);
530 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
556 sigma += mode2k2*bone2;
557 sigma += e2dk3t*bzero*bone;
560 sigma += kr2*bonebyarg2;
578 theta = std::sqrt(alpha);
582 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
595 bzero2 = bzero*bzero;
599 bonebyarg2 = bonebyarg*bonebyarg;
603 diffuse = 0.63*
fermi;
612 diffuse = 0.63*
fermi;
623 diffuse = 0.63*
fermi;
638 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
663 sigma += mode2k2*bone2;
664 sigma += e2dk3t*bzero*bone;
667 sigma += kr2*bonebyarg2;
724 G4double t = 2*p*p*( 1 - std::cos(theta) );
748 if (thetaMax >
pi) thetaMax =
pi;
757 for(i = 1; i <= iMax; i++)
759 theta1 = (i-1)*thetaMax/iMax;
760 theta2 = i*thetaMax/iMax;
765 result = 0.5*(theta1 + theta2);
769 if (i > iMax ) result = 0.5*(theta1 + theta2);
775 if(result < 0.) result = 0.;
776 if(result > thetaMax) result = thetaMax;
792 G4double totElab = std::sqrt(m1*m1+p*p);
807 G4double pCMS2 = momentumCMS*momentumCMS;
808 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
843 G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );
858 G4int iMomentum, iAngle;
880 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882 for( iMomentum = 0; iMomentum <
fEnergyBin; iMomentum++)
884 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
886 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
887 if ( iMomentum < 0 ) iMomentum = 0;
891 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
897 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
899 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
901 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
916 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
919 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
921 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
938 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
941 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
943 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
957 randAngle = W1*theta1 + W2*theta2;
965 if(randAngle < 0.) randAngle = 0.;
983 G4cout<<
"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
984 <<Z<<
"; and A = "<<A<<
G4endl;
1013 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1024 alphaMax = kRmax*kRmax/kR2;
1036 alphaCoulomb = kRcoul*kRcoul/kR2;
1041 fBeta = a/std::sqrt(1+a*a);
1065 alpha1 = delth*(j-1);
1067 alpha2 = alpha1 + delth;
1070 if( ( alpha1 < alphaCoulomb ) &&
z )
fAddCoulomb =
false;
1099 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1106 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1108 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1109 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1111 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1112 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1114 if ( x1 == x2 ) randAngle =
x2;
1120 randAngle = x1 + ( position -
y1 )*( x2 - x1 )/( y2 -
y1 );
1159 t =
SampleT( theParticle, ptot, A);
1162 if(!(t < 0.0 || t >= 0.0))
1166 G4cout <<
"G4DiffuseElastic:WARNING: A = " << A
1167 <<
" mom(GeV)= " << plab/
GeV
1168 <<
" S-wave will be sampled"
1175 G4cout <<
" t= " << t <<
" tmax= " << tmax
1176 <<
" ptot= " << ptot <<
G4endl;
1189 else if( cost <= -1.0)
1196 sint = std::sqrt((1.0-cost)*(1.0+cost));
1200 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1243 G4double cost = std::cos(thetaCMS);
1251 else if( cost <= -1.0)
1258 sint = std::sqrt((1.0-cost)*(1.0+cost));
1262 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1303 G4double cost = std::cos(thetaLab);
1311 else if( cost <= -1.0)
1318 sint = std::sqrt((1.0-cost)*(1.0+cost));
1322 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1352 G4cout<<
"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1353 <<Z<<
"; and A = "<<A<<
G4endl;
1362 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1363 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1364 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1378 alphaMax = kRmax*kRmax/kR2;
1380 if (alphaMax > 4.) alphaMax = 4.;
1382 alphaCoulomb = kRcoul*kRcoul/kR2;
1387 fBeta = a/std::sqrt(1+a*a);
1403 alpha1 = alphaMax*(j-1)/fAngleBin;
1404 alpha2 = alphaMax*( j )/fAngleBin;
1406 if( ( alpha2 > alphaCoulomb ) &&
z )
fAddCoulomb =
true;
1411 alpha1, alpha2,epsilon);
1421 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1423 angleVector->
PutValue( j-1 , alpha1, sumL10 );