G4RKPropagation.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 //
00027 // -------------------------------------------------------------------
00028 //      GEANT 4 class implementation file
00029 //
00030 //      CERN, Geneva, Switzerland
00031 //
00032 //      File name:     G4RKPropagation.cc
00033 //
00034 //      Author:        Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
00035 //
00036 //      Creation date: 6 June 2000
00037 // -------------------------------------------------------------------
00038 #include "G4RKPropagation.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 // nuclear fields
00042 #include "G4VNuclearField.hh"
00043 #include "G4ProtonField.hh"
00044 #include "G4NeutronField.hh"
00045 #include "G4AntiProtonField.hh"
00046 #include "G4KaonPlusField.hh"
00047 #include "G4KaonMinusField.hh"
00048 #include "G4KaonZeroField.hh"
00049 #include "G4PionPlusField.hh"
00050 #include "G4PionMinusField.hh"
00051 #include "G4PionZeroField.hh"
00052 #include "G4SigmaPlusField.hh"
00053 #include "G4SigmaMinusField.hh"
00054 #include "G4SigmaZeroField.hh"
00055 // particles properties
00056 #include "G4Proton.hh"
00057 #include "G4Neutron.hh"
00058 #include "G4AntiProton.hh"
00059 #include "G4KaonPlus.hh"
00060 #include "G4KaonMinus.hh"
00061 #include "G4KaonZero.hh"
00062 #include "G4PionPlus.hh"
00063 #include "G4PionMinus.hh"
00064 #include "G4PionZero.hh"
00065 #include "G4SigmaPlus.hh"
00066 #include "G4SigmaMinus.hh"
00067 #include "G4SigmaZero.hh"
00068 
00069 #include "globals.hh"
00070 
00071 #include "G4KM_OpticalEqRhs.hh"
00072 #include "G4KM_NucleonEqRhs.hh"
00073 #include "G4ClassicalRK4.hh"
00074 #include "G4MagIntegratorDriver.hh"
00075 
00076 #include "G4LorentzRotation.hh"
00077 
00078 // unsigned EncodingHashFun(const G4int& aEncoding);
00079 
00080 G4RKPropagation::G4RKPropagation() :
00081 theOuterRadius(0), theNucleus(0),
00082 theFieldMap(0), theEquationMap(0),
00083 theField(0)
00084 { }
00085 
00086 
00087 G4RKPropagation::~G4RKPropagation()
00088 {
00089    // free theFieldMap memory
00090    if(theFieldMap) delete_FieldsAndMap(theFieldMap);
00091 
00092    // free theEquationMap memory
00093    if(theEquationMap) delete_EquationsAndMap(theEquationMap);
00094 
00095    if (theField) delete theField;
00096 }
00097 
00098 //----------------------------------------------------------------------------
00099 void G4RKPropagation::Init(G4V3DNucleus * nucleus)
00100 //----------------------------------------------------------------------------
00101 {
00102 
00103    // free theFieldMap memory
00104    if(theFieldMap) delete_FieldsAndMap(theFieldMap);
00105 
00106    // free theEquationMap memory
00107    if(theEquationMap) delete_EquationsAndMap(theEquationMap);
00108 
00109    if (theField) delete theField;
00110 
00111    // Initialize the nuclear field map.
00112    theNucleus = nucleus;
00113    theOuterRadius = theNucleus->GetOuterRadius();
00114 
00115    theFieldMap = new std::map <G4int, G4VNuclearField*, std::less<G4int> >;
00116 
00117    (*theFieldMap)[G4Proton::Proton()->GetPDGEncoding()] = new G4ProtonField(theNucleus);
00118    (*theFieldMap)[G4Neutron::Neutron()->GetPDGEncoding()] = new G4NeutronField(theNucleus);
00119    (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = new G4AntiProtonField(theNucleus);
00120    (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = new G4KaonPlusField(theNucleus);
00121    (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = new G4KaonMinusField(theNucleus);
00122    (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = new G4KaonZeroField(theNucleus);
00123    (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = new G4PionPlusField(theNucleus);
00124    (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = new G4PionMinusField(theNucleus);
00125    (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()] = new G4PionZeroField(theNucleus);
00126    (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = new G4SigmaPlusField(theNucleus);
00127    (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = new G4SigmaMinusField(theNucleus);
00128    (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = new G4SigmaZeroField(theNucleus);
00129 
00130    theEquationMap = new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >;
00131 
00132    // theField needed by the design of G4Mag_eqRhs
00133    theField = new G4KM_DummyField;              //Field not needed for integration
00134    G4KM_OpticalEqRhs * opticalEq;
00135    G4KM_NucleonEqRhs * nucleonEq;
00136    G4double mass;
00137    G4double opticalCoeff;
00138 
00139    nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
00140    mass = G4Proton::Proton()->GetPDGMass();
00141    nucleonEq->SetMass(mass);
00142    (*theEquationMap)[G4Proton::Proton()->GetPDGEncoding()] = nucleonEq;
00143 
00144    nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
00145    mass = G4Neutron::Neutron()->GetPDGMass();
00146    nucleonEq->SetMass(mass);
00147    (*theEquationMap)[G4Neutron::Neutron()->GetPDGEncoding()] = nucleonEq;
00148 
00149    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00150    mass = G4AntiProton::AntiProton()->GetPDGMass();
00151    opticalCoeff =
00152          (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()]->GetCoeff();
00153    opticalEq->SetFactor(mass,opticalCoeff);
00154    (*theEquationMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = opticalEq;
00155 
00156    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00157    mass = G4KaonPlus::KaonPlus()->GetPDGMass();
00158    opticalCoeff =
00159          (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()]->GetCoeff();
00160    opticalEq->SetFactor(mass,opticalCoeff);
00161    (*theEquationMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = opticalEq;
00162 
00163    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00164    mass = G4KaonMinus::KaonMinus()->GetPDGMass();
00165    opticalCoeff =
00166          (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()]->GetCoeff();
00167    opticalEq->SetFactor(mass,opticalCoeff);
00168    (*theEquationMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = opticalEq;
00169 
00170    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00171    mass = G4KaonZero::KaonZero()->GetPDGMass();
00172    opticalCoeff =
00173          (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()]->GetCoeff();
00174    opticalEq->SetFactor(mass,opticalCoeff);
00175    (*theEquationMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = opticalEq;
00176 
00177    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00178    mass = G4PionPlus::PionPlus()->GetPDGMass();
00179    opticalCoeff =
00180          (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()]->GetCoeff();
00181    opticalEq->SetFactor(mass,opticalCoeff);
00182    (*theEquationMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = opticalEq;
00183 
00184    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00185    mass = G4PionMinus::PionMinus()->GetPDGMass();
00186    opticalCoeff =
00187          (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()]->GetCoeff();
00188    opticalEq->SetFactor(mass,opticalCoeff);
00189    (*theEquationMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = opticalEq;
00190 
00191    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00192    mass = G4PionZero::PionZero()->GetPDGMass();
00193    opticalCoeff =
00194          (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()]->GetCoeff();
00195    opticalEq->SetFactor(mass,opticalCoeff);
00196    (*theEquationMap)[G4PionZero::PionZero()->GetPDGEncoding()] = opticalEq;
00197 
00198    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00199    mass = G4SigmaPlus::SigmaPlus()->GetPDGMass();
00200    opticalCoeff =
00201          (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()]->GetCoeff();
00202    opticalEq->SetFactor(mass,opticalCoeff);
00203    (*theEquationMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = opticalEq;
00204 
00205    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00206    mass = G4SigmaMinus::SigmaMinus()->GetPDGMass();
00207    opticalCoeff =
00208          (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()]->GetCoeff();
00209    opticalEq->SetFactor(mass,opticalCoeff);
00210    (*theEquationMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = opticalEq;
00211 
00212    opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
00213    mass = G4SigmaZero::SigmaZero()->GetPDGMass();
00214    opticalCoeff =
00215          (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()]->GetCoeff();
00216    opticalEq->SetFactor(mass,opticalCoeff);
00217    (*theEquationMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = opticalEq;
00218 }
00219 
00220 
00221 //#define debug_1_RKPropagation 1
00222 //----------------------------------------------------------------------------
00223 void G4RKPropagation::Transport(G4KineticTrackVector & active,
00224       //----------------------------------------------------------------------------
00225       const G4KineticTrackVector &,
00226       G4double timeStep)
00227 {
00228    //  reset momentum transfer to field
00229    theMomentumTranfer=G4ThreeVector(0,0,0);
00230 
00231    // Loop over tracks
00232 
00233    std::vector<G4KineticTrack *>::iterator i;
00234    for(i = active.begin(); i != active.end(); ++i)
00235    {
00236       G4double currTimeStep = timeStep;
00237       G4KineticTrack * kt = *i;
00238       G4int encoding = kt->GetDefinition()->GetPDGEncoding();
00239 
00240       std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
00241 
00242       G4VNuclearField* currentField=0;
00243       if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second;
00244 
00245       // debug
00246       //    if ( timeStep > 1e30 ) {
00247       //        G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl;
00248       //    }
00249 
00250       // Get the time of intersections with the nucleus surface.
00251       G4double t_enter, t_leave;
00252       // if the particle does not intersecate with the nucleus go to next particle
00253       if(!GetSphereIntersectionTimes(kt, t_enter, t_leave))
00254       {
00255          kt->SetState(G4KineticTrack::miss_nucleus);
00256          continue;
00257       }
00258 
00259 
00260 #ifdef debug_1_RKPropagation
00261       G4cout <<" kt,timeStep, Intersection times tenter, tleave  "
00262             <<kt<< " / state= " <<kt->GetState() <<" / " <<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl;
00263 #endif
00264 
00265       // if the particle is already outside nucleus go to next  @@GF should never happen? check!
00266       //  does happen for particles added as late....
00267       //     if(t_leave < 0 )
00268       //     {
00269       //        throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a  nucleus");
00270       //        continue;
00271       //     }
00272 
00273       // Apply a straight line propagation for particle types
00274       // not included in the model
00275       if( ! currentField )
00276       {
00277          if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05;
00278          FreeTransport(kt, currTimeStep);
00279          if ( currTimeStep >= t_leave )
00280          {
00281             if ( kt->GetState() == G4KineticTrack::inside )
00282             { kt->SetState(G4KineticTrack::gone_out); }
00283             else
00284             { kt->SetState(G4KineticTrack::miss_nucleus);}
00285          } else if (kt->GetState() == G4KineticTrack::outside && currTimeStep >= t_enter ){
00286             kt->SetState(G4KineticTrack::inside);
00287          }
00288 
00289          continue;
00290       }
00291 
00292       if(t_enter > 0)  // the particle is out. Transport free to the surface
00293       {
00294          if(t_enter > currTimeStep)  // the particle won't enter the nucleus
00295          {
00296             FreeTransport(kt, currTimeStep);
00297             continue;
00298          }
00299          else
00300          {
00301             FreeTransport(kt, t_enter);   // go to surface
00302             currTimeStep -= t_enter;
00303             t_leave          -= t_enter;  // time left to leave nucleus
00304             // on the surface the particle loose the barrier energy
00305             //  G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier();
00306             //     GetField = Barrier + FermiPotential
00307             G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition());
00308 
00309             if(newE <= kt->GetActualMass())  // the particle cannot enter the nucleus
00310             {
00311                // FixMe: should be "pushed back?"
00312                //      for the moment take it past the nucleus, so we'll not worry next time..
00313                FreeTransport(kt, 1.1*t_leave);   // take past nucleus
00314                kt->SetState(G4KineticTrack::miss_nucleus);
00315                //          G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl;
00316                //          G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
00317                //          G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
00318                //          G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
00319                continue;
00320             }
00321             //
00322             G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
00323             G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
00324             G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
00325             G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
00326             new4Mom*=G4LorentzRotation(boost);
00327             kt->SetTrackingMomentum(new4Mom);
00328             kt->SetState(G4KineticTrack::inside);
00329 
00330             /*
00331      G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
00332            << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
00333            << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
00334            << G4endl
00335            << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
00336            << (*theFieldMap)[encoding]->GetBarrier() << " / "
00337            << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition())
00338            << G4endl;
00339              */
00340          }
00341       }
00342 
00343       // FixMe: should I add a control on theCutOnP here?
00344       // Transport the particle into the nucleus
00345       //       G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl;
00346       G4bool is_exiting=false;
00347       if(currTimeStep > t_leave)  // particle will exit from the nucleus
00348       {
00349          currTimeStep = t_leave;
00350          is_exiting=true;
00351       }
00352 
00353 #ifdef debug_1_RKPropagation
00354       G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl;
00355       G4cout << "RKPropagation Ekin, field, projectile potential, p "
00356             << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
00357             << kt->GetPosition()<<" "
00358             << G4endl << currentField->GetField(kt->GetPosition()) << " "
00359             <<  kt->GetProjectilePotential()<< G4endl
00360             << kt->GetTrackingMomentum()
00361             << G4endl;
00362 #endif
00363 
00364       G4LorentzVector momold=kt->GetTrackingMomentum();
00365       G4ThreeVector posold=kt->GetPosition();
00366 
00367       //    if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() ||
00368       if (currTimeStep > 0 &&
00369             ! FieldTransport(kt, currTimeStep)) {
00370          FreeTransport(kt,currTimeStep);
00371       }
00372 
00373 #ifdef debug_1_RKPropagation
00374       G4cout << "RKPropagation Ekin, field, p "
00375             << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
00376             << G4endl << currentField->GetField(kt->GetPosition())<< G4endl
00377             << kt->GetTrackingMomentum()
00378             << G4endl
00379             << "delta p " << momold-kt->GetTrackingMomentum() << G4endl
00380             << "del pos " << posold-kt->GetPosition()
00381             << G4endl;
00382 #endif
00383 
00384       // complete the transport
00385       // FixMe: in some cases there could be a significant
00386       //        part to do still in the nucleus, or we stepped to far... depending on
00387       //        slope of potential
00388       G4double t_in=-1, t_out=0;  // set onto boundary.
00389 
00390       // should go out, or are already out by a too long step..
00391       if(is_exiting ||
00392             (GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 ))  // particle is exiting
00393       {
00394          if(t_in < 0 && t_out >= 0)   //still inside, transport safely out.
00395          {
00396             // transport free to a position that is surely out of the nucleus, to avoid
00397             // a new transportation and a new adding the barrier next loop.
00398             G4ThreeVector savePos = kt->GetPosition();
00399             FreeTransport(kt, t_out);
00400             // and evaluate the right the energy
00401             G4double newE=kt->GetTrackingMomentum().e();
00402 
00403             //  G4cout << " V pos/savePos << "
00404             //          << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / "
00405             //          << (*theFieldMap)[encoding]->GetField(savePos)
00406             //          << G4endl;
00407 
00408             if ( std::abs(currentField->GetField(savePos)) > 0. &&
00409                   std::abs(currentField->GetField(kt->GetPosition())) > 0.)
00410             { // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0
00411                //           This wrongly adds or subtracts the Barrier here while
00412                //           this is done later.
00413                newE += currentField->GetField(savePos)
00414                                       - currentField->GetField(kt->GetPosition());
00415             }
00416 
00417             //       G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
00418 
00419             if(newE < kt->GetActualMass())
00420             {
00421 #ifdef debug_1_RKPropagation
00422                G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl;
00423                G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
00424 #endif
00425                if (kt->GetDefinition() == G4Proton::Proton() ||
00426                      kt->GetDefinition() == G4Neutron::Neutron() ) {
00427                   kt->SetState(G4KineticTrack::captured);
00428                } else {
00429                   kt->SetState(G4KineticTrack::gone_out);  //@@GF tofix
00430                }
00431                continue; // the particle cannot exit the nucleus
00432             }
00433             G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
00434             G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
00435             G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
00436             G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
00437             new4Mom*=G4LorentzRotation(boost);
00438             kt->SetTrackingMomentum(new4Mom);
00439          }
00440          // add the potential barrier
00441          // FixMe the Coulomb field is not parallel to mom, this is simple approximation
00442          G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition());
00443          if(newE < kt->GetActualMass())
00444          {  // the particle cannot exit the nucleus  @@@ GF check.
00445 #ifdef debug_1_RKPropagation
00446             G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
00447 #endif
00448             if (kt->GetDefinition() == G4Proton::Proton() ||
00449                   kt->GetDefinition() == G4Neutron::Neutron() ) {
00450                kt->SetState(G4KineticTrack::captured);
00451             } else {
00452                kt->SetState(G4KineticTrack::gone_out);  //@@GF tofix
00453             }
00454             continue;
00455          }
00456          G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
00457          G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
00458          G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
00459          G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
00460          new4Mom*=G4LorentzRotation(boost);
00461          kt->SetTrackingMomentum(new4Mom);
00462          kt->SetState(G4KineticTrack::gone_out);
00463       }
00464 
00465    }
00466 
00467 }
00468 
00469 
00470 //----------------------------------------------------------------------------
00471 G4ThreeVector G4RKPropagation::GetMomentumTransfer() const
00472 //----------------------------------------------------------------------------
00473 {
00474    return theMomentumTranfer;
00475 }
00476 
00477 
00478 //----------------------------------------------------------------------------
00479 G4bool G4RKPropagation::FieldTransport(G4KineticTrack * kt, const G4double timeStep)
00480 //----------------------------------------------------------------------------
00481 {
00482    theMomentumTranfer=G4ThreeVector(0,0,0);
00483    //    G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
00484    // create the integrator stepper
00485    //    G4Mag_EqRhs * equation = mapIter->second;
00486    G4Mag_EqRhs * equation = (*theEquationMap)[kt->GetDefinition()->GetPDGEncoding()];
00487    G4MagIntegratorStepper * stepper = new G4ClassicalRK4(equation);
00488 
00489    // create the integrator driver
00490    G4double hMin = 1.0e-25*second;   // arbitrary choice. Means 0.03 fm at c
00491    G4MagInt_Driver * driver = new G4MagInt_Driver(hMin, stepper);
00492 
00493    // Temporary: use driver->AccurateAdvance()
00494    // create the G4FieldTrack needed by AccurateAdvance
00495    G4double curveLength = 0;
00496    G4FieldTrack track(kt->GetPosition(),
00497          kt->GetTrackingMomentum().vect().unit(), // momentum direction
00498          curveLength, // curvelength
00499          kt->GetTrackingMomentum().e()-kt->GetActualMass(), // kinetic energy
00500          kt->GetActualMass(), // restmass
00501          kt->GetTrackingMomentum().beta()*c_light); // velocity
00502    // integrate
00503    G4double eps = 0.01;
00504    //    G4cout << "currTimeStep = " << currTimeStep << G4endl;
00505    if(!driver->AccurateAdvance(track, timeStep, eps))
00506    {  // cannot track this particle
00507 #ifdef debug_1_RKPropagation
00508       std::cerr << "G4RKPropagation::FieldTransport() warning: integration error."
00509             << G4endl << "position " << kt->GetPosition() << " 4mom " <<kt->GetTrackingMomentum()
00510             <<G4endl << " timestep " <<timeStep
00511             << G4endl;
00512 #endif
00513       delete driver;
00514       delete stepper;
00515       return false;
00516    }
00517    /*
00518        G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
00519            << (*theFieldMap)[encoding]->GetField(pos) << " / "
00520            << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
00521            << G4endl;
00522     */
00523 
00524    // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nuclues frame.
00525    G4ThreeVector MomentumTranfer = kt->GetTrackingMomentum().vect() - track.GetMomentum();
00526    G4ThreeVector boost= MomentumTranfer / std::sqrt (MomentumTranfer.mag2() +sqr(theNucleus->GetMass()));
00527 
00528    // update the kt
00529    kt->SetPosition(track.GetPosition());
00530    G4LorentzVector mom(track.GetMomentum(),std::sqrt(track.GetMomentum().mag2() + sqr(kt->GetActualMass())));
00531    mom *= G4LorentzRotation( boost );
00532    theMomentumTranfer += ( kt->GetTrackingMomentum() - mom ).vect();
00533    kt->SetTrackingMomentum(mom);
00534 
00535    //    G4cout <<"Stepper output"<<kt<<" "<<kt->GetTrackingMomentum()<<" "<<kt->GetPosition()<<G4endl;
00536    /*
00537     *   G4ThreeVector MomentumTranfer2=kt->GetTrackingMomentum().vect() - mom.vect();
00538     * G4cout << " MomentumTransfer/corrected" <<    MomentumTranfer << " " <<  MomentumTranfer.mag()
00539     *       <<  " " <<    MomentumTranfer2 << " " <<  MomentumTranfer2.mag() << " "
00540     *       << MomentumTranfer-MomentumTranfer2 << " "<<
00541     *       MomentumTranfer-MomentumTranfer2.mag() << " " << G4endl;
00542     *      G4cout <<" E/Field/Sum aft : " <<mom.e() << " / "
00543     *            << " / " << (*theFieldMap)[encoding]->GetField(pos)<< " / "
00544     *      << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
00545     *      << G4endl;
00546     */
00547 
00548    delete driver;
00549    delete stepper;
00550    return true;
00551 }
00552 
00553 //----------------------------------------------------------------------------
00554 G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep)
00555 //----------------------------------------------------------------------------
00556 {
00557    G4ThreeVector newpos = kt->GetPosition() +
00558          timeStep*c_light/kt->GetTrackingMomentum().e() * kt->GetTrackingMomentum().vect();
00559    kt->SetPosition(newpos);
00560    return true;
00561 }
00562 
00563 /*
00564 G4bool G4RKPropagation::WillBeCaptured(const G4KineticTrack * kt)
00565 {
00566   G4double radius = theOuterRadius;
00567 
00568 // evaluate the final energy. Il will be captured if newE or newP < 0
00569   G4ParticleDefinition * definition = kt->GetDefinition();
00570   G4double mass = definition->GetPDGMass();
00571   G4ThreeVector pos = kt->GetPosition();
00572   G4LorentzVector mom = kt->GetTrackingMomentum();
00573   G4VNuclearField * field = (*theFieldMap)[definition->GetPDGEncoding()];
00574   G4ThreeVector newPos(0, 0, radius); // to get the field on the surface
00575 
00576   G4double newE = mom.e()+field->GetField(pos)-field->GetField(newPos);
00577 
00578   return ((newE < mass) ? false : true);
00579 }
00580  */
00581 
00582 
00583 
00584 //----------------------------------------------------------------------------
00585 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4double radius,
00586       //----------------------------------------------------------------------------
00587       const G4ThreeVector & currentPos,
00588       const G4LorentzVector & momentum,
00589       G4double & t1, G4double & t2)
00590 {
00591    G4ThreeVector speed = momentum.vect()/momentum.e(); // boost vector
00592    G4double scalarProd = currentPos.dot(speed);
00593    G4double speedMag2 = speed.mag2();
00594    G4double sqrtArg = scalarProd*scalarProd -
00595          speedMag2*(currentPos.mag2()-radius*radius);
00596    if(sqrtArg <= 0.) // particle will not intersect the sphere
00597    {
00598       //     G4cout << " GetSphereIntersectionTimes sqrtArg negative: " << sqrtArg << G4endl;
00599       return false;
00600    }
00601    t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
00602    t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
00603    return true;
00604 }
00605 
00606 //----------------------------------------------------------------------------
00607 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4KineticTrack * kt,
00608       G4double & t1, G4double & t2)
00609 {
00610    G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi
00611    G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector
00612    G4double scalarProd = kt->GetPosition().dot(speed);
00613    G4double speedMag2 = speed.mag2();
00614    G4double sqrtArg = scalarProd*scalarProd -
00615          speedMag2*(kt->GetPosition().mag2()-radius*radius);
00616    if(sqrtArg <= 0.) // particle will not intersect the sphere
00617    {
00618       return false;
00619    }
00620    t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
00621    t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
00622    return true;
00623 }
00624 
00625 // Implementation methods
00626 
00627 //----------------------------------------------------------------------------
00628 void G4RKPropagation::delete_FieldsAndMap(
00629       //----------------------------------------------------------------------------
00630       std::map <G4int, G4VNuclearField *, std::less<G4int> > * aMap)
00631 {
00632    if(aMap)
00633    {
00634       std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur;
00635       for(cur = aMap->begin(); cur != aMap->end(); ++cur)
00636          delete (*cur).second;
00637 
00638       aMap->clear();
00639       delete aMap;
00640    }
00641 
00642 }
00643 
00644 //----------------------------------------------------------------------------
00645 void G4RKPropagation::delete_EquationsAndMap(
00646       //----------------------------------------------------------------------------
00647       std::map <G4int, G4Mag_EqRhs *, std::less<G4int> > * aMap)
00648 {
00649    if(aMap)
00650    {
00651       std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur;
00652       for(cur = aMap->begin(); cur != aMap->end(); ++cur)
00653          delete (*cur).second;
00654 
00655       aMap->clear();
00656       delete aMap;
00657    }
00658 }

Generated on Mon May 27 17:49:44 2013 for Geant4 by  doxygen 1.4.7