G4BinaryLightIonReaction Class Reference

#include <G4BinaryLightIonReaction.hh>

Inheritance diagram for G4BinaryLightIonReaction:

G4HadronicInteraction

Public Member Functions

 G4BinaryLightIonReaction (G4VPreCompoundModel *ptr=0)
virtual ~G4BinaryLightIonReaction ()
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
void SetPrecompound (G4VPreCompoundModel *ptr)
void SetDeExcitation (G4ExcitationHandler *ptr)
virtual void ModelDescription (std::ostream &) const

Detailed Description

Definition at line 34 of file G4BinaryLightIonReaction.hh.


Constructor & Destructor Documentation

G4BinaryLightIonReaction::G4BinaryLightIonReaction ( G4VPreCompoundModel ptr = 0  ) 

Definition at line 49 of file G4BinaryLightIonReaction.cc.

References G4HadronicInteractionRegistry::FindModel(), G4VPreCompoundModel::GetExcitationHandler(), and G4HadronicInteractionRegistry::Instance().

00050 : G4HadronicInteraction("Binary Light Ion Cascade"),
00051   theProjectileFragmentation(ptr),
00052   pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
00053   projectile3dNucleus(0),target3dNucleus(0)
00054 {
00055         if(!ptr) {
00056           G4HadronicInteraction* p =
00057             G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
00058           G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
00059           if(!pre) { pre = new G4PreCompoundModel(); }
00060           theProjectileFragmentation = pre;
00061         }
00062         theModel = new G4BinaryCascade(theProjectileFragmentation);
00063         theHandler = theProjectileFragmentation->GetExcitationHandler();
00064 
00065         debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0;
00066 }

G4BinaryLightIonReaction::~G4BinaryLightIonReaction (  )  [virtual]

Definition at line 68 of file G4BinaryLightIonReaction.cc.

00069 {}


Member Function Documentation

G4HadFinalState * G4BinaryLightIonReaction::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus theNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 87 of file G4BinaryLightIonReaction.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cerr, G4cout, G4endl, G4lrint(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4HadFinalState::GetNumberOfSecondaries(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4HadProjectile::GetTotalEnergy(), G4Nucleus::GetZ_asInt(), isAlive, G4DynamicParticle::Set4Momentum(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), and stopAndKill.

00088 {
00089         static G4int eventcounter=0;
00090         eventcounter++;
00091         if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number starts ######### "<<eventcounter<<G4endl;
00092         G4ping debug("debug_G4BinaryLightIonReaction");
00093         pA=aTrack.GetDefinition()->GetBaryonNumber();
00094         pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus);
00095         tA=targetNucleus.GetA_asInt();
00096         tZ=targetNucleus.GetZ_asInt();
00097 
00098         G4LorentzVector mom(aTrack.Get4Momentum());
00099    //G4cout << "proj mom : " << mom << G4endl;
00100         G4LorentzRotation toBreit(mom.boostVector());
00101 
00102         G4bool swapped=SetLighterAsProjectile(mom, toBreit);
00103    //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
00104         G4ReactionProductVector * result = 0;
00105         G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
00106 //      G4double m_nucl(0);      // to check energy balance
00107 
00108 
00109         //    G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
00110         //    G4cout << "Entering the decision point "
00111         //           << (mom.t()-mom.mag())/pA << " "
00112         //         << pA<<" "<< pZ<<" "
00113         //         << tA<<" "<< tZ<<G4endl
00114         //         << " "<<mom.t()-mom.mag()<<" "
00115         //         << mom.t()- m1<<G4endl;
00116         if( (mom.t()-mom.mag())/pA < 50*MeV )
00117         {
00118                 //      G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
00119                 //      m_nucl = mom.mag();
00120       cascaders=FuseNucleiAndPrompound(mom);
00121         }
00122         else
00123         {
00124            result=Interact(mom,toBreit);
00125 
00126       if(! result )
00127       {
00128           {
00129             // abort!!
00130 
00131             G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
00132             G4cerr << " Primary " << aTrack.GetDefinition()
00133                << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
00134                << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
00135                << ", kinetic energy " << aTrack.GetKineticEnergy()
00136                << G4endl;
00137             G4cerr << " Target nucleus (A,Z)=("
00138                    <<  (swapped?pA:tA)  << ","
00139                    << (swapped?pZ:tZ) << ")" << G4endl;
00140             G4cerr << " if frequent, please submit above information as bug report"
00141                   << G4endl << G4endl;
00142 
00143             theResult.Clear();
00144             theResult.SetStatusChange(isAlive);
00145             theResult.SetEnergyChange(aTrack.GetKineticEnergy());
00146             theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00147             return &theResult;
00148 
00149          }
00150       }
00151 
00152                 // Calculate excitation energy,
00153       G4double theStatisticalExEnergy = GetProjectileExcitation();
00154 
00155 
00156                 pInitialState = mom;
00157                 //G4cout << "pInitialState from aTrack : " << pInitialState;
00158                 pInitialState.setT(pInitialState.getT() +
00159                      G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA));
00160       //G4cout << "target nucleus added : " << pInitialState << G4endl;
00161 
00162                 delete target3dNucleus;target3dNucleus=0;
00163                 delete projectile3dNucleus;projectile3dNucleus=0;
00164 
00165       G4ReactionProductVector * spectators= new G4ReactionProductVector;
00166 
00167       cascaders = new G4ReactionProductVector;
00168 
00169       G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
00170 
00171       //      pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
00172       std::vector<G4ReactionProduct *>::iterator iter;
00173 
00174       //G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
00175                 //      if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
00176                 //      {
00177                 //          G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
00178                 //            spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
00179                 //      }
00180                 delete result;
00181                 result=0;
00182                 G4LorentzVector momentum(pInitialState-pFinalState);
00183                 G4int loopcount(0);
00184                 //G4cout << "momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
00185                 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV)
00186                 {
00187                         G4LorentzVector pCorrect(pInitialState-pspectators);
00188                    //G4cout << "BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
00189                         // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
00190                         G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
00191                         if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
00192                         {
00193                                 G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
00194                         }
00195                         pFinalState=G4LorentzVector(0,0,0,0);
00196                         unsigned int i;
00197                         for(i=0; i<cascaders->size(); i++)
00198                         {
00199                                 pFinalState += G4LorentzVector( (*cascaders)[i]->GetMomentum(), (*cascaders)[i]->GetTotalEnergy() );
00200                         }
00201                         momentum=pInitialState-pFinalState;
00202                         if (++loopcount > 10 )
00203                         {
00204                                 if ( momentum.vect().mag() - momentum.e()> 10*keV  )
00205                                 {
00206                                         G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl;
00207                                         throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
00208                                 } else {
00209                                         break;
00210                                 }
00211 
00212                         }
00213                 }
00214 
00215         if (spectatorA > 0 )
00216                 {
00217                    // check spectator momentum
00218                    if ( momentum.vect().mag() - momentum.e()> 10*keV )
00219                    {
00220 
00221                       G4ReactionProductVector::iterator ispectator;
00222                       for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
00223                       {
00224                          delete *ispectator;
00225                       }
00226                       delete spectators;
00227 
00228                       G4cout << "G4BinaryLightIonReaction.cc: mom check: " <<  momentum
00229                             << " 3.mag "<< momentum.vect().mag() << G4endl
00230                             << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
00231                             << pFinalState << " " << pspectators << G4endl
00232                             << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
00233                       G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
00234                       G4cout << " Primary " << aTrack.GetDefinition()
00235                               << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
00236                               << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
00237                               << ", kinetic energy " << aTrack.GetKineticEnergy()
00238                               << G4endl;
00239                       G4cout << " Target nucleus (A,Z)=(" <<  targetNucleus.GetA_asInt()
00240                               << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
00241                       G4cout << " if frequent, please submit above information as bug report"
00242                             << G4endl << G4endl;
00243 
00244                       theResult.Clear();
00245                       theResult.SetStatusChange(isAlive);
00246                       theResult.SetEnergyChange(aTrack.GetKineticEnergy());
00247                       theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00248                       return &theResult;
00249                    }
00250 
00251 
00252                    DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
00253                 }
00254         }
00255         // Rotate to lab
00256         G4LorentzRotation toZ;
00257         toZ.rotateZ(-1*mom.phi());
00258         toZ.rotateY(-1*mom.theta());
00259         G4LorentzRotation toLab(toZ.inverse());
00260 
00261         // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
00262         // theResult.Clear();
00263         theResult.Clear();
00264         theResult.SetStatusChange(stopAndKill);
00265         G4double Etot(0);
00266         size_t i=0;
00267         for(i=0; i<cascaders->size(); i++)
00268         {
00269                 if((*cascaders)[i]->GetNewlyAdded())
00270                 {
00271                         G4DynamicParticle * aNew =
00272                                         new G4DynamicParticle((*cascaders)[i]->GetDefinition(),
00273                                                         (*cascaders)[i]->GetTotalEnergy(),
00274                                                         (*cascaders)[i]->GetMomentum() );
00275                         G4LorentzVector tmp = aNew->Get4Momentum();
00276                         if(swapped)
00277                         {
00278                                 tmp*=toBreit.inverse();
00279                                 tmp.setVect(-tmp.vect());
00280                         }
00281                         tmp *= toLab;
00282                         aNew->Set4Momentum(tmp);
00283                         //G4cout << "result[" << i << "], 4vect: " << tmp << G4endl;
00284                         theResult.AddSecondary(aNew);
00285                         Etot += tmp.e();
00286                         //        G4cout << "LIBIC: Secondary " << aNew->GetDefinition()->GetParticleName()
00287                         //               <<" "<<  aNew->GetMomentum()
00288                         //            <<" "<<  aNew->GetTotalEnergy()
00289                         //            << G4endl;
00290                 }
00291                 delete (*cascaders)[i];
00292         }
00293         delete cascaders;
00294 
00295 #ifdef debug_BLIR_result
00296         G4cout << "Result analysis, secondaries" << theResult.GetNumberOfSecondaries() << G4endl;
00297         G4cout << " Energy conservation initial/primary/nucleus/final/delta(init-final) "
00298               << aTrack.GetTotalEnergy() + m_nucl << aTrack.GetTotalEnergy() << m_nucl <<Etot
00299               << aTrack.GetTotalEnergy() + m_nucl - Etot;
00300 #endif
00301 
00302         if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### "<<eventcounter<<G4endl;
00303 
00304         return &theResult;
00305 }

void G4BinaryLightIonReaction::ModelDescription ( std::ostream &   )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 71 of file G4BinaryLightIonReaction.cc.

00072 {
00073         outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
00074                         << "using G4BinaryCasacde to model the interaction of a light\n"
00075                         << "nucleus with a nucleus.\n"
00076                         << "The lighter of the two nuclei is treated like a set of projectiles\n"
00077                         << "which are transported simultanously through the heavier nucleus.\n";
00078 }

void G4BinaryLightIonReaction::SetDeExcitation ( G4ExcitationHandler ptr  )  [inline]

Definition at line 74 of file G4BinaryLightIonReaction.hh.

References G4VPreCompoundModel::SetExcitationHandler().

00075 {
00076   theProjectileFragmentation->SetExcitationHandler(ptr);
00077   theHandler = ptr;
00078 }

void G4BinaryLightIonReaction::SetPrecompound ( G4VPreCompoundModel ptr  )  [inline]

Definition at line 69 of file G4BinaryLightIonReaction.hh.

References G4VPreCompoundModel::GetExcitationHandler().

00070 {
00071   if(ptr) { theProjectileFragmentation = ptr; }
00072   theHandler = theProjectileFragmentation->GetExcitationHandler();
00073 }


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