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 #include <vector>
00030
00031 #include "globals.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4ios.hh"
00035 #include "G4Scatterer.hh"
00036 #include "G4KineticTrack.hh"
00037 #include "G4ThreeVector.hh"
00038 #include "G4LorentzRotation.hh"
00039 #include "G4LorentzVector.hh"
00040
00041 #include "G4CollisionNN.hh"
00042 #include "G4CollisionPN.hh"
00043 #include "G4CollisionMesonBaryon.hh"
00044
00045 #include "G4CollisionInitialState.hh"
00046 #include "G4HadTmpUtil.hh"
00047 #include "G4Pair.hh"
00048
00049
00050 typedef GROUP2(G4CollisionNN, G4CollisionMesonBaryon) theChannels;
00051
00052
00053
00054 G4Scatterer::G4Scatterer()
00055 {
00056 Register aR;
00057 G4ForEach<theChannels>::Apply(&aR, &collisions);
00058 }
00059
00060
00061
00062 G4Scatterer::~G4Scatterer()
00063 {
00064 std::for_each(collisions.begin(), collisions.end(), G4Delete());
00065 collisions.clear();
00066 }
00067
00068
00069
00070 G4double G4Scatterer::GetTimeToInteraction(const G4KineticTrack& trk1,
00071 const G4KineticTrack& trk2)
00072 {
00073 G4double time = DBL_MAX;
00074 G4double distance_fast;
00075 G4LorentzVector mom1 = trk1.GetTrackingMomentum();
00076
00077 G4double collisionTime;
00078
00079 if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
00080 {
00081 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
00082 G4double deltaz=position.z();
00083 G4double velocity = mom1.z()/mom1.e() * c_light;
00084
00085 collisionTime=deltaz/velocity;
00086 distance_fast=position.x()*position.x() + position.y()*position.y();
00087 } else {
00088
00089
00090
00091 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
00092
00093 G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;
00094 collisionTime = (position * velocity) / velocity.mag2();
00095 position -= velocity * collisionTime;
00096 distance_fast=position.mag2();
00097
00098
00099
00100 }
00101 if (collisionTime > 0)
00102 {
00103 static const G4double maxCrossSection = 500*millibarn;
00104 if(0.7*pi*distance_fast>maxCrossSection) return time;
00105
00106
00107 G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
00108
00109
00110
00111
00112
00113 G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
00114 mom1 = toCMSFrame * mom1;
00115 mom2 = toCMSFrame * mom2;
00116
00117 G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
00118 G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
00119 G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
00120 (toCMSFrame * coordinate2).vect());
00121
00122 G4ThreeVector mom = mom1.vect() - mom2.vect();
00123
00124
00125
00126 G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
00127
00128
00129
00130
00131
00132
00133
00134 if(pi*distance>maxCrossSection) return time;
00135
00136
00137 static const G4double maxChargedCrossSection = 200*millibarn;
00138 if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
00139 std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
00140 pi*distance>maxChargedCrossSection) return time;
00141
00142 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00143
00144 if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
00145 trk2.GetDefinition() == G4Neutron::Neutron() ) &&
00146 sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181 if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
00182
00183
00184
00185 G4VCollision* collision = FindCollision(trk1,trk2);
00186 G4double totalCrossSection;
00187
00188
00189
00190 if (collision != 0)
00191 {
00192 totalCrossSection = collision->CrossSection(trk1,trk2);
00193 if ( totalCrossSection > 0 )
00194 {
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 if (distance <= totalCrossSection / pi)
00207 {
00208 time = collisionTime;
00209 }
00210 } else
00211 {
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 }
00236 else
00237 {
00238 time = DBL_MAX;
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 }
00249 }
00250
00251 else
00252 {
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 }
00263
00264 return time;
00265 }
00266
00267
00268
00269 G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1,
00270 const G4KineticTrack& trk2)
00271 {
00272
00273 G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
00274 G4double energyBalance = pInitial.t();
00275 G4double pxBalance = pInitial.vect().x();
00276 G4double pyBalance = pInitial.vect().y();
00277 G4double pzBalance = pInitial.vect().z();
00278 G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
00279 + trk2.GetDefinition()->GetPDGCharge());
00280 G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
00281 + trk2.GetDefinition()->GetBaryonNumber();
00282
00283 G4VCollision* collision = FindCollision(trk1,trk2);
00284 if (collision != 0)
00285 {
00286 G4double aCrossSection = collision->CrossSection(trk1,trk2);
00287 if (aCrossSection > 0.0)
00288 {
00289
00290
00291 #ifdef debug_G4Scatterer
00292 G4cout << "be4 FinalState 1(p,e,m): "
00293 << trk1.Get4Momentum() << " "
00294 << trk1.Get4Momentum().mag()
00295 << ", 2: "
00296 << trk2.Get4Momentum()<< " "
00297 << trk2.Get4Momentum().mag() << " "
00298 << G4endl;
00299 #endif
00300
00301
00302 G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
00303 if(!products || products->size() == 0) return products;
00304
00305 #ifdef debug_G4Scatterer
00306 G4cout << "size of FS: "<<products->size()<<G4endl;
00307 #endif
00308
00309 G4KineticTrack *final= products->operator[](0);
00310
00311
00312 #ifdef debug_G4Scatterer
00313 G4cout << " FinalState 1: "
00314 << final->Get4Momentum()<< " "
00315 << final->Get4Momentum().mag() ;
00316 #endif
00317
00318 if(products->size() == 1) return products;
00319 final=products->operator[](1);
00320 #ifdef debug_G4Scatterer
00321 G4cout << ", 2: "
00322 << final->Get4Momentum() << " "
00323 << final->Get4Momentum().mag() << " " << G4endl;
00324 #endif
00325
00326 final= products->operator[](0);
00327 G4LorentzVector pFinal=final->Get4Momentum();
00328 if(products->size()==2)
00329 {
00330 final=products->operator[](1);
00331 pFinal +=final->Get4Momentum();
00332 }
00333
00334 #ifdef debug_G4Scatterer
00335 if ( (pInitial-pFinal).mag() > 0.1*MeV )
00336 {
00337 G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
00338 }
00339 G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
00340 #endif
00341
00342 for(size_t hpw=0; hpw<products->size(); hpw++)
00343 {
00344 energyBalance-=products->operator[](hpw)->Get4Momentum().t();
00345 pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
00346 pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
00347 pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
00348 chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
00349 baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
00350 }
00351 if(getenv("ScattererEnergyBalanceCheck"))
00352 std::cout << "DEBUGGING energy balance A: "
00353 <<energyBalance<<" "
00354 <<pxBalance<<" "
00355 <<pyBalance<<" "
00356 <<pzBalance<<" "
00357 <<chargeBalance<<" "
00358 <<baryonBalance<<" "
00359 <<G4endl;
00360 if(chargeBalance !=0 )
00361 {
00362 G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
00363 G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
00364 for(size_t hpw=0; hpw<products->size(); hpw++)
00365 {
00366 G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
00367 }
00368 G4Exception("G4Scatterer", "im_r_matrix001", FatalException,
00369 "Problem in ChargeBalance");
00370 }
00371 return products;
00372 }
00373 }
00374
00375 return NULL;
00376 }
00377
00378
00379
00380 G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1,
00381 const G4KineticTrack& trk2)
00382 {
00383 G4VCollision* collisionInCharge = 0;
00384
00385 size_t i;
00386 for (i=0; i<collisions.size(); i++)
00387 {
00388 G4VCollision* component = collisions[i];
00389 if (component->IsInCharge(trk1,trk2))
00390 {
00391 collisionInCharge = component;
00392 break;
00393 }
00394 }
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 return collisionInCharge;
00405 }
00406
00407
00408
00409 G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1,
00410 const G4KineticTrack& trk2)
00411 {
00412 G4VCollision* collision = FindCollision(trk1,trk2);
00413 G4double aCrossSection = 0;
00414 if (collision != 0)
00415 {
00416 aCrossSection = collision->CrossSection(trk1,trk2);
00417 }
00418 return aCrossSection;
00419 }
00420
00421
00422
00423 const std::vector<G4CollisionInitialState *> & G4Scatterer::
00424 GetCollisions(G4KineticTrack * aProjectile,
00425 std::vector<G4KineticTrack *> & someCandidates,
00426 G4double aCurrentTime)
00427 {
00428 theCollisions.clear();
00429 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
00430 for(; j != someCandidates.end(); ++j)
00431 {
00432 G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
00433 if(collisionTime == DBL_MAX)
00434 {
00435 continue;
00436 }
00437 G4KineticTrackVector aTarget;
00438 aTarget.push_back(*j);
00439 theCollisions.push_back(
00440 new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
00441
00442 }
00443 return theCollisions;
00444 }
00445
00446
00447 G4KineticTrackVector * G4Scatterer::
00448 GetFinalState(G4KineticTrack * aProjectile,
00449 std::vector<G4KineticTrack *> & theTargets)
00450 {
00451 G4KineticTrack target_reloc(*(theTargets[0]));
00452 return Scatter(*aProjectile, target_reloc);
00453 }
00454