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 #include "globals.hh"
00040 #include "Randomize.hh"
00041 #include "G4PhysicalConstants.hh"
00042
00043 #include "G4ElasticHNScattering.hh"
00044 #include "G4LorentzRotation.hh"
00045 #include "G4ThreeVector.hh"
00046 #include "G4ParticleDefinition.hh"
00047 #include "G4VSplitableHadron.hh"
00048 #include "G4ExcitedString.hh"
00049 #include "G4FTFParameters.hh"
00050
00051
00052 G4ElasticHNScattering::G4ElasticHNScattering()
00053 {
00054 }
00055
00056 G4bool G4ElasticHNScattering::
00057 ElasticScattering (G4VSplitableHadron *projectile,
00058 G4VSplitableHadron *target,
00059 G4FTFParameters *theParameters) const
00060 {
00061
00062 G4LorentzVector Pprojectile=projectile->Get4Momentum();
00063
00064 if(Pprojectile.z() < 0.)
00065 {
00066 target->SetStatus(2);
00067 return false;
00068 }
00069
00070 G4bool PutOnMassShell(false);
00071
00072 G4double M0projectile = Pprojectile.mag();
00073 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
00074 {
00075 PutOnMassShell=true;
00076 M0projectile=projectile->GetDefinition()->GetPDGMass();
00077 }
00078
00079 G4double Mprojectile2 = M0projectile * M0projectile;
00080
00081 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
00082
00083
00084
00085 G4LorentzVector Ptarget=target->Get4Momentum();
00086
00087 G4double M0target = Ptarget.mag();
00088
00089 if(M0target < target->GetDefinition()->GetPDGMass())
00090 {
00091 PutOnMassShell=true;
00092 M0target=target->GetDefinition()->GetPDGMass();
00093 }
00094
00095 G4double Mtarget2 = M0target * M0target;
00096
00097
00098
00099 G4LorentzVector Psum;
00100 Psum=Pprojectile+Ptarget;
00101
00102 G4LorentzRotation toCms(-1*Psum.boostVector());
00103
00104 G4LorentzVector Ptmp=toCms*Pprojectile;
00105
00106 if ( Ptmp.pz() <= 0. )
00107 {
00108
00109
00110 target->SetStatus(2);
00111 return false;
00112 }
00113
00114 toCms.rotateZ(-1*Ptmp.phi());
00115 toCms.rotateY(-1*Ptmp.theta());
00116
00117 G4LorentzRotation toLab(toCms.inverse());
00118
00119 Pprojectile.transform(toCms);
00120 Ptarget.transform(toCms);
00121
00122
00123 G4double PZcms2, PZcms;
00124
00125 G4double S=Psum.mag2();
00126
00127
00128 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00129 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
00130
00131 if(PZcms2 < 0.)
00132 {
00133 if(M0projectile > projectile->GetDefinition()->GetPDGMass())
00134 {
00135
00136 M0projectile = projectile->GetDefinition()->GetPDGMass();
00137 Mprojectile2=M0projectile*M0projectile;
00138 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00139 2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
00140 /4./S;
00141
00142 if(PZcms2 < 0.){ return false;}
00143 }
00144 else
00145 {
00146 target->SetStatus(2);
00147 return false;
00148
00149
00150 }
00151 }
00152
00153 PZcms = std::sqrt(PZcms2);
00154
00155 if(PutOnMassShell)
00156 {
00157 if(Pprojectile.z() > 0.)
00158 {
00159 Pprojectile.setPz( PZcms);
00160 Ptarget.setPz( -PZcms);
00161 }
00162 else
00163 {
00164 Pprojectile.setPz(-PZcms);
00165 Ptarget.setPz( PZcms);
00166 };
00167
00168 Pprojectile.setE(std::sqrt(Mprojectile2+
00169 Pprojectile.x()*Pprojectile.x()+
00170 Pprojectile.y()*Pprojectile.y()+
00171 PZcms2));
00172 Ptarget.setE(std::sqrt( Mtarget2 +
00173 Ptarget.x()*Ptarget.x()+
00174 Ptarget.y()*Ptarget.y()+
00175 PZcms2));
00176 }
00177
00178 G4double maxPtSquare = PZcms2;
00179
00180
00181 G4double Pt2;
00182 G4double ProjMassT2;
00183 G4double TargMassT2;
00184
00185 G4LorentzVector Qmomentum;
00186 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00187
00188 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00189
00190 ProjMassT2=Mprojectile2+Pt2;
00191
00192
00193 TargMassT2=Mtarget2+Pt2;
00194
00195
00196 PZcms2=(S*S+ProjMassT2*ProjMassT2+
00197 TargMassT2*TargMassT2-
00198 2.*S*ProjMassT2-2.*S*TargMassT2-
00199 2.*ProjMassT2*TargMassT2)/4./S;
00200 if(PZcms2 < 0 ) {PZcms2=0;};
00201 PZcms =std::sqrt(PZcms2);
00202
00203 Pprojectile.setPz( PZcms);
00204 Ptarget.setPz( -PZcms);
00205
00206 Pprojectile += Qmomentum;
00207 Ptarget -= Qmomentum;
00208
00209
00210 Pprojectile.transform(toLab);
00211 Ptarget.transform(toLab);
00212
00213
00214
00215
00216
00217
00218
00219 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00220 projectile->SetPosition(target->GetPosition());
00221
00222
00223
00224
00225 projectile->Set4Momentum(Pprojectile);
00226 target->Set4Momentum(Ptarget);
00227
00228 projectile->IncrementCollisionCount(1);
00229 target->IncrementCollisionCount(1);
00230
00231 return true;
00232 }
00233
00234
00235
00236
00237 G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
00238 {
00239
00240 G4double Pt2(0.);
00241 if(AveragePt2 <= 0.) {Pt2=0.;}
00242 else
00243 {
00244 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
00245 (std::exp(-maxPtSquare/AveragePt2)-1.));
00246 }
00247 G4double Pt=std::sqrt(Pt2);
00248 G4double phi=G4UniformRand() * twopi;
00249
00250 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
00251 }
00252
00253 G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
00254 {
00255 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
00256 }
00257
00258
00259 G4ElasticHNScattering::~G4ElasticHNScattering()
00260 {
00261 }
00262
00263
00264 const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
00265 {
00266 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator not meant to be called");
00267
00268 }
00269
00270
00271 int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
00272 {
00273 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator not meant to be called");
00274
00275 }
00276
00277 int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
00278 {
00279 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator not meant to be called");
00280
00281 }