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 #include "G4Absorber.hh"
00027 #include "G4KineticTrack.hh"
00028 #include "G4PionPlus.hh"
00029 #include "G4PionMinus.hh"
00030 #include "G4PionZero.hh"
00031 #include "G4Proton.hh"
00032 #include "G4Neutron.hh"
00033
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4LorentzRotation.hh"
00037
00038 G4Absorber::G4Absorber(G4double cutOnP)
00039 {
00040 theCutOnP = cutOnP;
00041 theAbsorbers = new G4KineticTrackVector;
00042 theProducts = new G4KineticTrackVector;
00043 }
00044
00045
00046 G4Absorber::~G4Absorber()
00047 {
00048 delete theAbsorbers;
00049 delete theProducts;
00050 }
00051
00052
00053 bool G4Absorber::WillBeAbsorbed(const G4KineticTrack & kt)
00054 {
00055
00056
00057
00058 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP)
00059 {
00060 if(kt.GetDefinition() == G4PionPlus::PionPlus() ||
00061 kt.GetDefinition() == G4PionZero::PionZero() ||
00062 kt.GetDefinition() == G4PionMinus::PionMinus())
00063 {
00064 return true;
00065 }
00066 }
00067 return false;
00068 }
00069
00070
00071
00072 G4bool G4Absorber::Absorb(G4KineticTrack & kt, G4KineticTrackVector & tgt)
00073 {
00074 if(!FindAbsorbers(kt, tgt))
00075 return false;
00076 return FindProducts(kt);
00077 }
00078
00079
00080 G4bool G4Absorber::FindAbsorbers(G4KineticTrack & kt,
00081 G4KineticTrackVector & tgt)
00082 {
00083
00084
00085
00086
00087
00088
00089 G4KineticTrack * kt1 = NULL;
00090 G4KineticTrack * kt2 = NULL;
00091 G4double dist1 = DBL_MAX;
00092 G4double dist2 = DBL_MAX;
00093 G4double charge1 = 0;
00094
00095 G4double charge0 = kt.GetDefinition()->GetPDGCharge();
00096 G4ThreeVector pos = kt.GetPosition();
00097
00098 std::vector<G4KineticTrack *>::iterator iter;
00099 for(iter = tgt.begin(); iter != tgt.end(); ++iter)
00100 {
00101 G4KineticTrack * curr = *iter;
00102 G4double dist = (pos-curr->GetPosition()).mag();
00103 if(dist >= dist2)
00104 continue;
00105 if(dist < dist1)
00106 {
00107 if(dist1 == DBL_MAX)
00108 {
00109 kt1 = curr;
00110 charge1 = kt1->GetDefinition()->GetPDGCharge();
00111 dist1 = dist;
00112 continue;
00113 }
00114 if(dist2 == DBL_MAX)
00115 {
00116 kt2 = kt1;
00117
00118 dist2 = dist1;
00119 kt1 = curr;
00120 charge1 = kt1->GetDefinition()->GetPDGCharge();
00121 dist1 = dist;
00122 continue;
00123 }
00124
00125 G4double charge = curr->GetDefinition()->GetPDGCharge();
00126 if((charge0+charge1+charge < 0.) ||
00127 (charge0+charge1+charge) > 2*eplus)
00128 {
00129 kt1 = curr;
00130 charge1 = charge;
00131 dist1 = dist;
00132 }
00133 else
00134 {
00135 kt2 = kt1;
00136
00137 dist2 = dist1;
00138 kt1 = curr;
00139 charge1 = charge;
00140 dist1 = dist;
00141 }
00142 continue;
00143 }
00144
00145 if(dist2 == DBL_MAX)
00146 {
00147 kt2 = curr;
00148
00149 dist2 = dist;
00150 continue;
00151 }
00152
00153 G4double charge = curr->GetDefinition()->GetPDGCharge();
00154 if((charge0+charge1+charge < 0.) ||
00155 (charge0+charge1+charge) > 2*eplus)
00156 continue;
00157
00158 kt2 = curr;
00159
00160 dist2 = dist;
00161 }
00162
00163 theAbsorbers->clear();
00164 if((kt1 == NULL) || (kt2 == NULL))
00165 return false;
00166
00167 theAbsorbers->push_back(kt1);
00168 theAbsorbers->push_back(kt2);
00169 return true;
00170 }
00171
00172
00173
00174 G4bool G4Absorber::FindProducts(G4KineticTrack & kt)
00175 {
00176
00177 G4ParticleDefinition * prod1;
00178 G4ParticleDefinition * prod2;
00179 G4KineticTrack * abs1 = (*theAbsorbers)[0];
00180 G4KineticTrack * abs2 = (*theAbsorbers)[1];
00181
00182 G4double charge = kt.GetDefinition()->GetPDGCharge();
00183 if(charge == eplus)
00184 {
00185 prod1 = G4Proton::Proton();
00186 if(abs1->GetDefinition() == G4Neutron::Neutron())
00187 prod2 = abs2->GetDefinition();
00188 else
00189 prod2 = G4Proton::Proton();
00190 }
00191 else if(charge == -eplus)
00192 {
00193 prod1 = G4Neutron::Neutron();
00194 if(abs1->GetDefinition() == G4Proton::Proton())
00195 prod2 = abs2->GetDefinition();
00196 else
00197 prod2 = G4Neutron::Neutron();
00198 }
00199 else
00200 {
00201 prod1 = abs1->GetDefinition();
00202 prod2 = abs2->GetDefinition();
00203 }
00204
00205
00206 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+
00207 abs2->Get4Momentum();
00208 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector());
00209 G4LorentzRotation toLabFrame(momLab.boostVector());
00210 G4LorentzVector momCMS = toCMSFrame*momLab;
00211
00212
00213 G4double ms1 = prod1->GetPDGMass();
00214 G4double ms2 = prod2->GetPDGMass();
00215 G4double e0 = momCMS.e();
00216 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1*ms1+ms2*ms2)+
00217 (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0*e0);
00218
00219
00220 G4ThreeVector mom1CMS = GetRandomDirection();
00221 mom1CMS = std::sqrt(squareP)*mom1CMS;
00222 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+ms1*ms1));
00223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+ms2*ms2));
00224
00225
00226 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS;
00227 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS;
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(),
00246 mom1);
00247 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(),
00248 mom2);
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 theProducts->clear();
00279 theProducts->push_back(kt1);
00280 theProducts->push_back(kt2);
00281 return true;
00282 }
00283
00284
00285
00286 G4ThreeVector G4Absorber::GetRandomDirection()
00287 {
00288 G4double theta = 2.0*G4UniformRand()-1.0;
00289 theta = std::acos(theta);
00290 G4double phi = G4UniformRand()*2*pi;
00291 G4ThreeVector direction(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
00292 return direction;
00293 }
00294
00295
00296
00297
00298
00299
00300