G4Scatterer Class Reference

#include <G4Scatterer.hh>

Inheritance diagram for G4Scatterer:

Public Member Functions
	G4Scatterer ()
virtual	~G4Scatterer ()
virtual G4double	GetTimeToInteraction (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
G4double	GetCrossSection (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
virtual G4KineticTrackVector *	Scatter (const G4KineticTrack &trk1, const G4KineticTrack &trk2)
virtual const std::vector < G4CollisionInitialState * > &	GetCollisions (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &someCandidates, G4double aCurrentTime)
virtual G4KineticTrackVector *	GetFinalState (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &theTargets)
Public Member Functions inherited from G4VScatterer
	G4VScatterer ()
virtual	~G4VScatterer ()
Public Member Functions inherited from G4BCAction
	G4BCAction ()
virtual	~G4BCAction ()

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(), and G4AutoDelete::Register().

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

G4Scatterer::~G4Scatterer ( )

virtual

Definition at line 62 of file G4Scatterer.cc.

{
  std::for_each(collisions.begin(), collisions.end(), G4Delete());
  collisions.clear();
}

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().

{
  theCollisions.clear();
  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
  for(; j != someCandidates.end(); ++j)
  {
    G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
    if(collisionTime == DBL_MAX)  // no collision
    {
      continue;
    }
    G4KineticTrackVector aTarget;
    aTarget.push_back(*j);
    theCollisions.push_back(
      new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
//      G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
   }
   return theCollisions;
}

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

Definition at line 409 of file G4Scatterer.cc.

References G4VCollision::CrossSection().

{
   G4VCollision* collision = FindCollision(trk1,trk2);
   G4double aCrossSection = 0;
   if (collision != 0)
   {
     aCrossSection = collision->CrossSection(trk1,trk2);
   }
   return aCrossSection;
}

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

virtual

Implements G4BCAction.

Definition at line 448 of file G4Scatterer.cc.

References Scatter().

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

G4double G4Scatterer::GetTimeToInteraction ( const G4KineticTrack &  trk1, const G4KineticTrack &  trk2  )

virtual

Implements G4VScatterer.

Definition at line 70 of file G4Scatterer.cc.

References python.hepunit::c_light, G4VCollision::CrossSection(), DBL_MAX, CLHEP::HepLorentzVector::e(), G4KineticTrack::Get4Momentum(), G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), G4KineticTrack::GetPosition(), G4KineticTrack::GetTrackingMomentum(), python.hepunit::GeV, CLHEP::HepLorentzVector::mag(), CLHEP::Hep3Vector::mag2(), python.hepunit::millibarn, G4Neutron::Neutron(), python.hepunit::pi, CLHEP::Hep3Vector::unit(), CLHEP::HepLorentzVector::vect(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), CLHEP::Hep3Vector::z(), and CLHEP::HepLorentzVector::z().

Referenced by GetCollisions().

{
  G4double time = DBL_MAX;
    G4double distance_fast;
  G4LorentzVector mom1 = trk1.GetTrackingMomentum();
//  G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
  G4double collisionTime;

  if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
  {
     G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
     G4double deltaz=position.z();
     G4double velocity = mom1.z()/mom1.e() * c_light;

     collisionTime=deltaz/velocity;
     distance_fast=position.x()*position.x() + position.y()*position.y();
  } else {

    //  The nucleons of the nucleus are FROZEN, ie. do not move..

    G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();

    G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;  // mom1.boostVector() will exit on slightly negative mass
    collisionTime = (position * velocity) / velocity.mag2();    // can't divide by /c_light;
    position -= velocity * collisionTime;
    distance_fast=position.mag2();

//    if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
//     collisionTime = GetTimeToClosestApproach(trk1,trk2);
  }
     if (collisionTime > 0)
    {
       static const G4double maxCrossSection = 500*millibarn;
       if(0.7*pi*distance_fast>maxCrossSection) return time;


           G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());

//     G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
//     G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
//     G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());

       G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
       mom1 = toCMSFrame * mom1;
       mom2 = toCMSFrame * mom2;

       G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
       G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
       G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
                (toCMSFrame * coordinate2).vect());

       G4ThreeVector mom = mom1.vect() - mom2.vect();

      // Calculate the impact parameter

       G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());

//     G4cout << " disDiff " << distance-disLab << " " << disLab
//            << " " << std::abs(distance-disLab)/distance << G4endl
//      << " mom/Lab " << mom << " " << momLab << G4endl
//      << " pos/Lab " << pos << " " << posLab
//      << G4endl;

       if(pi*distance>maxCrossSection) return time;

       // charged particles special
       static const G4double maxChargedCrossSection = 200*millibarn;
       if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
          std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
          pi*distance>maxChargedCrossSection) return time;

           G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
       // neutrons special  pn is largest cross-section, but above 1.91 GeV is less than 200 mb
       if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
            trk2.GetDefinition() == G4Neutron::Neutron() ) &&
             sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;

/*
 *    if(distance <= sqr(1.14*fermi))
 *    {
 *      time = collisionTime;
 *
 * *
 *  *        G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
 *  *            " / "<< time/ns << G4endl;
 *  *         G4ThreeVector pos1=trk1.GetPosition();
 *  *         G4ThreeVector pos2=trk2.GetPosition();
 *  *         G4LorentzVector xmom1 = trk1.Get4Momentum();
 *  *         G4LorentzVector xmom2 = trk2.Get4Momentum();
 *  *         G4cout << "position1: " <<  pos1.x() << " " << pos1.y() << " "
 *  *               << pos1.z();
 *  *         pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
 *  *         G4cout << " straight line trprt: "
 *  *               <<  pos1.x() << " " << pos1.y() << " "
 *  *           <<  pos1.z()  << G4endl;
 *  *         G4cout << "position2: " <<  pos2.x() << " " << pos2.y() << " "
 *  *               << pos2.z()  << G4endl;
 *  *         G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
 *  *         pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
 *  *         G4cout<< " straight line trprt: "
 *  *               <<  pos2.x() << " " << pos2.y() << " "
 *  *           <<  pos2.z() << G4endl;
 *  *         G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
 *  *
 *    }
 *
 *    if(1)
 *      return time;
 */

       if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;



      G4VCollision* collision = FindCollision(trk1,trk2);
      G4double totalCrossSection;
      // The cross section is interpreted geometrically as an area
      // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)

      if (collision != 0)
        {
          totalCrossSection = collision->CrossSection(trk1,trk2);
          if ( totalCrossSection > 0 )
            {
/*          G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
 *                 << trk1.GetDefinition()->GetParticleName()
 *             << " / "
 *                 << trk2.GetDefinition()->GetParticleName()
 *             << ", "
 *             << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
 *             << ", "
 *             << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
 *                 trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
 *             << G4endl;
 */
         if (distance <= totalCrossSection / pi)
           {
             time = collisionTime;
           }
            } else
        {

         // For debugging...
 //         G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
 //                << trk1.GetDefinition()->GetParticleName()
 //            << " / "
 //                << trk2.GetDefinition()->GetParticleName()
 //            << ", "
 //            << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
 //            << ", "
 //            << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
 //                trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
 //            << G4endl;

        }
/*
 *        if(distance <= sqr(5.*fermi))
 *         {
 *        G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
 *           << " " << totalCrossSection/sqr(fermi)
 *           << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
 *         }
 */

        }
      else
        {
          time = DBL_MAX;
//        /*
          // For debugging
//hpw                 G4cout << "G4Scatterer - collision not found: "
//hpw            << trk1.GetDefinition()->GetParticleName()
//hpw            << " - "
//hpw            << trk2.GetDefinition()->GetParticleName()
//hpw            << G4endl;
          // End of debugging
//        */
        }
    }

      else
    {
          /*
      // For debugging
      G4cout << "G4Scatterer - negative collisionTime"
         << ": collisionTime = " << collisionTime
         << ", position = " << position
         << ", velocity = " << velocity
         << G4endl;
      // End of debugging
          */
    }

  return time;
}

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(), G4ParticleDefinition::GetPDGCharge(), CLHEP::HepLorentzVector::mag(), python.hepunit::MeV, CLHEP::HepLorentzVector::t(), CLHEP::Hep3Vector::unit(), CLHEP::HepLorentzVector::vect(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

{
//   G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
   G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
   G4double energyBalance = pInitial.t();
   G4double pxBalance = pInitial.vect().x();
   G4double pyBalance = pInitial.vect().y();
   G4double pzBalance = pInitial.vect().z();
   G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
                       + trk2.GetDefinition()->GetPDGCharge());
   G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
                       + trk2.GetDefinition()->GetBaryonNumber();

   G4VCollision* collision = FindCollision(trk1,trk2);
   if (collision != 0)
   {
     G4double aCrossSection = collision->CrossSection(trk1,trk2);
     if (aCrossSection > 0.0)
     {


    #ifdef debug_G4Scatterer
    G4cout << "be4 FinalState 1(p,e,m): "
        << trk1.Get4Momentum() << " "
        << trk1.Get4Momentum().mag()
    << ", 2: "
        << trk2.Get4Momentum()<< " "
        << trk2.Get4Momentum().mag() << " "
        << G4endl;
    #endif


       G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
       if(!products || products->size() == 0) return products;

    #ifdef debug_G4Scatterer
       G4cout << "size of FS: "<<products->size()<<G4endl;
    #endif

       G4KineticTrack *final= products->operator[](0);


    #ifdef debug_G4Scatterer
        G4cout << "    FinalState 1: "
        << final->Get4Momentum()<< " "
        << final->Get4Momentum().mag() ;
    #endif

        if(products->size() == 1) return products;
    final=products->operator[](1);
    #ifdef debug_G4Scatterer
    G4cout << ", 2: "
        << final->Get4Momentum() << " "
            << final->Get4Momentum().mag() << " " << G4endl;
    #endif

       final= products->operator[](0);
       G4LorentzVector pFinal=final->Get4Momentum();
       if(products->size()==2)
       {
         final=products->operator[](1);
         pFinal +=final->Get4Momentum();
       }

       #ifdef debug_G4Scatterer
       if ( (pInitial-pFinal).mag() > 0.1*MeV )
       {
          G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
       }
       G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
       #endif

       for(size_t hpw=0; hpw<products->size(); hpw++)
       {
         energyBalance-=products->operator[](hpw)->Get4Momentum().t();
         pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
         pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
         pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
     chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
         baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
       }
       if(getenv("ScattererEnergyBalanceCheck"))
         std::cout << "DEBUGGING energy balance A: "
               <<energyBalance<<" "
               <<pxBalance<<" "
               <<pyBalance<<" "
               <<pzBalance<<" "
           <<chargeBalance<<" "
           <<baryonBalance<<" "
           <<G4endl;
       if(chargeBalance !=0 )
       {
         G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
         G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
         for(size_t hpw=0; hpw<products->size(); hpw++)
         {
            G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
         }
         G4Exception("G4Scatterer", "im_r_matrix001", FatalException,
             "Problem in ChargeBalance");
       }
       return products;
     }
   }

   return NULL;
}

---

The documentation for this class was generated from the following files:

• G4Scatterer.hh
• G4Scatterer.cc