G4Absorber Class Reference

#include <G4Absorber.hh>


Public Member Functions

 G4Absorber (G4double cutOnP)
 ~G4Absorber ()
G4bool WillBeAbsorbed (const G4KineticTrack &kt)
G4bool Absorb (G4KineticTrack &kt, G4KineticTrackVector &tgt)
G4KineticTrackVectorGetAbsorbers ()
G4KineticTrackVectorGetProducts ()
G4bool FindAbsorbers (G4KineticTrack &kt, G4KineticTrackVector &tgt)
G4bool FindProducts (G4KineticTrack &kt)


Detailed Description

Definition at line 34 of file G4Absorber.hh.


Constructor & Destructor Documentation

G4Absorber::G4Absorber ( G4double  cutOnP  ) 

Definition at line 38 of file G4Absorber.cc.

00039 {
00040   theCutOnP = cutOnP;
00041   theAbsorbers = new G4KineticTrackVector;
00042   theProducts = new G4KineticTrackVector;
00043 }

G4Absorber::~G4Absorber (  ) 

Definition at line 46 of file G4Absorber.cc.

00047 {
00048   delete theAbsorbers;
00049   delete theProducts;
00050 }


Member Function Documentation

G4bool G4Absorber::Absorb ( G4KineticTrack kt,
G4KineticTrackVector tgt 
)

Definition at line 72 of file G4Absorber.cc.

References FindAbsorbers(), and FindProducts().

00073 {
00074   if(!FindAbsorbers(kt, tgt))
00075     return false;
00076   return FindProducts(kt);
00077 }

G4bool G4Absorber::FindAbsorbers ( G4KineticTrack kt,
G4KineticTrackVector tgt 
)

Definition at line 80 of file G4Absorber.cc.

References DBL_MAX, G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), and G4KineticTrack::GetPosition().

Referenced by Absorb().

00082 {
00083 //  Find a closest ( in space) pair of Nucleons capable to absorb pi+/pi-
00084 //    pi+ can be absorbed on np or nn resulting in pp or np
00085 //    pi- can be absorbed on np or pp resulting in nn or np
00086 
00087 // @GF: FindAbsorbers is unused, logic is seriously wrong
00088 
00089   G4KineticTrack * kt1 = NULL;
00090   G4KineticTrack * kt2 = NULL;
00091   G4double dist1 = DBL_MAX;             // dist to closest nucleon
00092   G4double dist2 = DBL_MAX;             // dist to next close 
00093   G4double charge1 = 0;
00094 //  G4double charge2 = 0;               // charge2 is only assigned to, never used
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) // accept 1st as a candidate, 
00108       {
00109         kt1 = curr;
00110         charge1 = kt1->GetDefinition()->GetPDGCharge();
00111         dist1 = dist;
00112         continue;
00113       } 
00114       if(dist2 == DBL_MAX) // accept the candidate and shift kt1 to kt2
00115       {                         // @GF: should'nt we check if compatible?
00116         kt2 = kt1;
00117 //      charge2 = charge1;
00118         dist2 = dist1;
00119         kt1 = curr;
00120         charge1 = kt1->GetDefinition()->GetPDGCharge();
00121         dist1 = dist;
00122         continue;
00123       }
00124 // test the compatibility with charge conservation for new config
00125       G4double charge = curr->GetDefinition()->GetPDGCharge();
00126       if((charge0+charge1+charge < 0.) ||       //test config (curr,kt1) 
00127          (charge0+charge1+charge) > 2*eplus)
00128       {  // incompatible: change kt1 with curr.
00129         kt1 = curr;
00130         charge1 = charge;
00131         dist1 = dist;
00132       }
00133       else
00134       { // compatible: change kt1 with curr and kt2 with kt1
00135         kt2 = kt1;
00136 //      charge2 = charge1;
00137         dist2 = dist1;
00138         kt1 = curr;
00139         charge1 = charge;
00140         dist1 = dist;
00141       }
00142       continue;
00143     }
00144 // here if dist1 < dist < dist2
00145     if(dist2 == DBL_MAX) // accept the candidate
00146     {
00147       kt2 = curr;
00148 //      charge2 = kt2->GetDefinition()->GetPDGCharge();
00149       dist2 = dist;
00150       continue;
00151     }   
00152 // test the compatibility with charge conservation
00153     G4double charge = curr->GetDefinition()->GetPDGCharge();
00154     if((charge0+charge1+charge < 0.) ||
00155        (charge0+charge1+charge) > 2*eplus)
00156       continue;   // incomatible: do nothing
00157 // compatible: change kt2 with curr
00158     kt2 = curr;
00159 //    charge2 = charge;
00160     dist2 = dist;
00161   }
00162 
00163   theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector!
00164   if((kt1 == NULL) || (kt2 == NULL))
00165     return false;
00166 
00167   theAbsorbers->push_back(kt1);
00168   theAbsorbers->push_back(kt2);
00169   return true;
00170 }

G4bool G4Absorber::FindProducts ( G4KineticTrack kt  ) 

Definition at line 174 of file G4Absorber.cc.

References G4KineticTrack::Get4Momentum(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4KineticTrack::GetPosition(), G4Neutron::Neutron(), and G4Proton::Proton().

Referenced by Absorb().

00175 {
00176 // Choose the products type
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   { // a neutron become proton
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   { // a proton become neutron
00193     prod1 = G4Neutron::Neutron();
00194     if(abs1->GetDefinition() == G4Proton::Proton())
00195       prod2 = abs2->GetDefinition();
00196     else
00197       prod2 = G4Neutron::Neutron();
00198   }
00199   else  // charge = 0: leave particle types unchenged
00200   {
00201     prod1 = abs1->GetDefinition();
00202     prod2 = abs2->GetDefinition();
00203   }
00204 
00205 // Translate to the CMS frame
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 // Evaluate the final momentum of products
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 //  if(squareP < 0)  // should never happen
00219 //    squareP = 0;
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 // Go back to the lab frame
00226   G4LorentzVector mom1 = toLabFrame*final4Mom1CMS;
00227   G4LorentzVector mom2 = toLabFrame*final4Mom2CMS;
00228 
00229 // ------ debug
00230 /*
00231   G4LorentzVector temp = mom1+mom2;
00232 
00233   cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " "
00234        << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " "
00235        << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " "
00236        << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " "
00237        << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " "
00238        << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " "
00239        << (1/MeV)*std::sqrt(squareP) << endl;
00240 
00241 */
00242 // ------ end debug
00243 
00244 // Build two new kinetic tracks and add to products
00245   G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(),
00246                                             mom1);
00247   G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(),
00248                                             mom2);
00249 // ------ debug
00250 /*
00251   G4LorentzVector initialMom1 = abs1->Get4Momentum();
00252   G4LorentzVector initialMom2 = abs2->Get4Momentum();
00253   G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum();
00254   cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " "
00255        << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " "
00256        << (1/MeV)*initialMom1.vect().mag() << " "
00257        << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " "
00258        << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " "
00259        << (1/MeV)*initialMom2.vect().mag() << " "
00260        << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " "
00261        << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " "
00262        << (1/MeV)*mom1.vect().mag() << " "
00263        << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " "
00264        << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " "
00265        << (1/MeV)*mom2.vect().mag() << " "
00266        << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " "
00267        << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " "
00268        << (1/MeV)*pion4MomCMS.vect().mag() << " "
00269        << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " "
00270        << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " "
00271        << (1/MeV)*final4Mom1CMS.vect().mag() << " "
00272        << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " "
00273        << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " "
00274        << (1/MeV)*final4Mom2CMS.vect().mag() << endl;
00275 */
00276 // ------ end debug
00277 
00278   theProducts->clear();
00279   theProducts->push_back(kt1);
00280   theProducts->push_back(kt2);
00281   return true;
00282 }

G4KineticTrackVector * G4Absorber::GetAbsorbers (  )  [inline]

Definition at line 66 of file G4Absorber.hh.

00067 {
00068   return theAbsorbers;
00069 }

G4KineticTrackVector * G4Absorber::GetProducts (  )  [inline]

Definition at line 71 of file G4Absorber.hh.

00072 {
00073   return theProducts;
00074 }

bool G4Absorber::WillBeAbsorbed ( const G4KineticTrack kt  ) 

Definition at line 53 of file G4Absorber.cc.

References G4KineticTrack::Get4Momentum(), G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), and G4PionZero::PionZero().

00054 {
00055  // FixMe: actually only for pions
00056 //  if(kt.Get4Momentum().vect().mag() < theCutOnP) 
00057 // Cut on kinetic Energy...
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 }


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