G4Absorber.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
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  // 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 }
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 //  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 }
00171 
00172 
00173 
00174 G4bool G4Absorber::FindProducts(G4KineticTrack & kt)
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 }
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 

Generated on Mon May 27 17:47:36 2013 for Geant4 by  doxygen 1.4.7