G4Scatterer Class Reference

#include <G4Scatterer.hh>

Inheritance diagram for G4Scatterer:

G4VScatterer G4BCAction

Public Member Functions

 G4Scatterer ()
virtual ~G4Scatterer ()
virtual G4double GetTimeToInteraction (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
G4double GetCrossSection (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
virtual G4KineticTrackVectorScatter (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
virtual const std::vector<
G4CollisionInitialState * > & 
GetCollisions (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &someCandidates, G4double aCurrentTime)
virtual G4KineticTrackVectorGetFinalState (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &theTargets)

Data Structures

struct  Register

Detailed Description

Definition at line 44 of file G4Scatterer.hh.


Constructor & Destructor Documentation

G4Scatterer::G4Scatterer (  ) 

Definition at line 54 of file G4Scatterer.cc.

References G4ForEach< group >::Apply().

00055 {
00056   Register aR;
00057   G4ForEach<theChannels>::Apply(&aR, &collisions);
00058 }

G4Scatterer::~G4Scatterer (  )  [virtual]

Definition at line 62 of file G4Scatterer.cc.

00063 {
00064   std::for_each(collisions.begin(), collisions.end(), G4Delete());
00065   collisions.clear();
00066 }


Member Function Documentation

const std::vector< G4CollisionInitialState * > & G4Scatterer::GetCollisions ( G4KineticTrack aProjectile,
std::vector< G4KineticTrack * > &  someCandidates,
G4double  aCurrentTime 
) [virtual]

Implements G4BCAction.

Definition at line 424 of file G4Scatterer.cc.

References DBL_MAX, and GetTimeToInteraction().

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)  // no collision
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 //      G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
00442    }
00443    return theCollisions;
00444 }

G4double G4Scatterer::GetCrossSection ( const G4KineticTrack trk1,
const G4KineticTrack trk2 
)

Definition at line 409 of file G4Scatterer.cc.

References G4VCollision::CrossSection().

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 }

G4KineticTrackVector * G4Scatterer::GetFinalState ( G4KineticTrack aProjectile,
std::vector< G4KineticTrack * > &  theTargets 
) [virtual]

Implements G4BCAction.

Definition at line 448 of file G4Scatterer.cc.

References Scatter().

00450 {
00451     G4KineticTrack target_reloc(*(theTargets[0]));
00452     return Scatter(*aProjectile, target_reloc);
00453 }

G4double G4Scatterer::GetTimeToInteraction ( const G4KineticTrack trk1,
const G4KineticTrack trk2 
) [virtual]

Implements G4VScatterer.

Definition at line 70 of file G4Scatterer.cc.

References G4VCollision::CrossSection(), DBL_MAX, G4KineticTrack::Get4Momentum(), G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), G4KineticTrack::GetPosition(), G4KineticTrack::GetTrackingMomentum(), G4Neutron::Neutron(), G4INCL::Math::pi, and position.

Referenced by GetCollisions().

00072 {
00073   G4double time = DBL_MAX;
00074     G4double distance_fast;
00075   G4LorentzVector mom1 = trk1.GetTrackingMomentum();
00076 //  G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
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     //  The nucleons of the nucleus are FROZEN, ie. do not move..
00090 
00091     G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
00092 
00093     G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;  // mom1.boostVector() will exit on slightly negative mass
00094     collisionTime = (position * velocity) / velocity.mag2();    // can't divide by /c_light;
00095     position -= velocity * collisionTime;
00096     distance_fast=position.mag2();
00097 
00098 //    if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
00099 //     collisionTime = GetTimeToClosestApproach(trk1,trk2);
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 //         G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
00110 //         G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
00111 //         G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
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           // Calculate the impact parameter
00125 
00126            G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
00127 
00128 //     G4cout << " disDiff " << distance-disLab << " " << disLab
00129 //            << " " << std::abs(distance-disLab)/distance << G4endl
00130 //          << " mom/Lab " << mom << " " << momLab << G4endl
00131 //          << " pos/Lab " << pos << " " << posLab
00132 //          << G4endl;
00133 
00134            if(pi*distance>maxCrossSection) return time;
00135 
00136            // charged particles special
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            // neutrons special  pn is largest cross-section, but above 1.91 GeV is less than 200 mb
00144            if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
00145                 trk2.GetDefinition() == G4Neutron::Neutron() ) &&
00146                      sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
00147 
00148 /*
00149  *        if(distance <= sqr(1.14*fermi))
00150  *        {
00151  *          time = collisionTime;
00152  *
00153  * *
00154  *  *        G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
00155  *  *            " / "<< time/ns << G4endl;
00156  *  *         G4ThreeVector pos1=trk1.GetPosition();
00157  *  *         G4ThreeVector pos2=trk2.GetPosition();
00158  *  *         G4LorentzVector xmom1 = trk1.Get4Momentum();
00159  *  *         G4LorentzVector xmom2 = trk2.Get4Momentum();
00160  *  *         G4cout << "position1: " <<  pos1.x() << " " << pos1.y() << " "
00161  *  *                   << pos1.z();
00162  *  *         pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
00163  *  *         G4cout << " straight line trprt: "
00164  *  *                   <<  pos1.x() << " " << pos1.y() << " "
00165  *  *                   <<  pos1.z()  << G4endl;
00166  *  *         G4cout << "position2: " <<  pos2.x() << " " << pos2.y() << " "
00167  *  *                   << pos2.z()  << G4endl;
00168  *  *         G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
00169  *  *         pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
00170  *  *         G4cout<< " straight line trprt: "
00171  *  *                   <<  pos2.x() << " " << pos2.y() << " "
00172  *  *                   <<  pos2.z() << G4endl;
00173  *  *         G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
00174  *  *
00175  *        }
00176  *
00177  *        if(1)
00178  *          return time;
00179  */
00180 
00181            if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
00182 
00183 
00184 
00185           G4VCollision* collision = FindCollision(trk1,trk2);
00186           G4double totalCrossSection;
00187           // The cross section is interpreted geometrically as an area
00188           // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
00189 
00190           if (collision != 0)
00191             {
00192               totalCrossSection = collision->CrossSection(trk1,trk2);
00193               if ( totalCrossSection > 0 )
00194                 {
00195 /*                  G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
00196  *                         << trk1.GetDefinition()->GetParticleName()
00197  *                         << " / "
00198  *                         << trk2.GetDefinition()->GetParticleName()
00199  *                         << ", "
00200  *                         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
00201  *                         << ", "
00202  *                         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
00203  *                             trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
00204  *                         << G4endl;
00205  */
00206                  if (distance <= totalCrossSection / pi)
00207                    {
00208                      time = collisionTime;
00209                    }
00210                 } else
00211                 {
00212 
00213                  // For debugging...
00214  //                 G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
00215  //                        << trk1.GetDefinition()->GetParticleName()
00216  //                        << " / "
00217  //                        << trk2.GetDefinition()->GetParticleName()
00218  //                        << ", "
00219  //                        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
00220  //                        << ", "
00221  //                        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
00222  //                            trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
00223  //                        << G4endl;
00224 
00225                 }
00226 /*
00227  *            if(distance <= sqr(5.*fermi))
00228  *             {
00229  *                G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
00230  *                       << " " << totalCrossSection/sqr(fermi)
00231  *                       << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
00232  *             }
00233  */
00234 
00235             }
00236           else
00237             {
00238               time = DBL_MAX;
00239 //            /*
00240               // For debugging
00241 //hpw                 G4cout << "G4Scatterer - collision not found: "
00242 //hpw                << trk1.GetDefinition()->GetParticleName()
00243 //hpw                << " - "
00244 //hpw                << trk2.GetDefinition()->GetParticleName()
00245 //hpw                << G4endl;
00246               // End of debugging
00247 //            */
00248             }
00249         }
00250 
00251       else
00252         {
00253               /*
00254           // For debugging
00255           G4cout << "G4Scatterer - negative collisionTime"
00256                  << ": collisionTime = " << collisionTime
00257                  << ", position = " << position
00258                  << ", velocity = " << velocity
00259                  << G4endl;
00260           // End of debugging
00261               */
00262         }
00263 
00264   return time;
00265 }

G4KineticTrackVector * G4Scatterer::Scatter ( const G4KineticTrack trk1,
const G4KineticTrack trk2 
) [virtual]

Implements G4VScatterer.

Definition at line 269 of file G4Scatterer.cc.

References G4VCollision::CrossSection(), FatalException, G4VCollision::FinalState(), G4cout, G4endl, G4Exception(), G4lrint(), G4KineticTrack::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetParticleName(), and G4ParticleDefinition::GetPDGCharge().

Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision(), and GetFinalState().

00271 {
00272 //   G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
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 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:20 2013 for Geant4 by  doxygen 1.4.7