G4NeutronHPInelasticCompFS Class Reference

#include <G4NeutronHPInelasticCompFS.hh>

Inheritance diagram for G4NeutronHPInelasticCompFS:

G4NeutronHPFinalState G4NeutronHPAInelasticFS G4NeutronHPDInelasticFS G4NeutronHPHe3InelasticFS G4NeutronHPNInelasticFS G4NeutronHPPInelasticFS G4NeutronHPTInelasticFS

Public Member Functions

 G4NeutronHPInelasticCompFS ()
virtual ~G4NeutronHPInelasticCompFS ()
void Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType)
void InitGammas (G4double AR, G4double ZR)
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &theTrack)=0
virtual G4NeutronHPFinalStateNew ()=0
virtual G4double GetXsec (G4double anEnergy)
virtual G4NeutronHPVectorGetXsec ()
G4int SelectExitChannel (G4double eKinetic)
void CompositeApply (const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)
void InitDistributionInitialState (G4ReactionProduct &aNeutron, G4ReactionProduct &aTarget, G4int it)

Protected Attributes

G4NeutronHPVectortheXsection [51]
G4NeutronHPEnergyDistributiontheEnergyDistribution [51]
G4NeutronHPAngulartheAngularDistribution [51]
G4NeutronHPEnAngCorrelationtheEnergyAngData [51]
G4NeutronHPPhotonDisttheFinalStatePhotons [51]
G4NeutronHPDeExGammas theGammas
G4String gammaPath
G4double theCurrentA
G4double theCurrentZ
std::vector< G4doubleQI
std::vector< G4intLR

Detailed Description

Definition at line 42 of file G4NeutronHPInelasticCompFS.hh.


Constructor & Destructor Documentation

G4NeutronHPInelasticCompFS::G4NeutronHPInelasticCompFS (  )  [inline]

Definition at line 46 of file G4NeutronHPInelasticCompFS.hh.

References G4NeutronHPFinalState::hasXsec, LR, QI, theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

00047   {
00048 
00049     QI.resize(51);
00050     LR.resize(51);
00051     for(G4int i=0; i<51; i++)
00052     {
00053       hasXsec = true; 
00054       theXsection[i] = 0;
00055       theEnergyDistribution[i] = 0;
00056       theAngularDistribution[i] = 0;
00057       theEnergyAngData[i] = 0;
00058       theFinalStatePhotons[i] = 0;
00059       QI[i]=0.0;
00060       LR[i]=0;
00061     }
00062 
00063   }

virtual G4NeutronHPInelasticCompFS::~G4NeutronHPInelasticCompFS (  )  [inline, virtual]

Definition at line 64 of file G4NeutronHPInelasticCompFS.hh.

References theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

00065   {
00066     for(G4int i=0; i<51; i++)
00067     {
00068       if(theXsection[i] != 0) delete theXsection[i];
00069       if(theEnergyDistribution[i] != 0) delete theEnergyDistribution[i];
00070       if(theAngularDistribution[i] != 0) delete theAngularDistribution[i];
00071       if(theEnergyAngData[i] != 0) delete theEnergyAngData[i];
00072       if(theFinalStatePhotons[i] != 0) delete theFinalStatePhotons[i];
00073     }
00074   }


Member Function Documentation

virtual G4HadFinalState* G4NeutronHPInelasticCompFS::ApplyYourself ( const G4HadProjectile theTrack  )  [pure virtual]

Reimplemented from G4NeutronHPFinalState.

Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

void G4NeutronHPInelasticCompFS::CompositeApply ( const G4HadProjectile theTrack,
G4ParticleDefinition aHadron 
)

Definition at line 216 of file G4NeutronHPInelasticCompFS.cc.

References G4HadFinalState::AddSecondary(), G4NeutronHPFinalState::adjust_final_state(), G4HadFinalState::Clear(), G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4Nucleus::GetBiasedThermalNucleus(), G4NeutronHPDeExGammas::GetDecayGammas(), G4ReactionProduct::GetDefinition(), G4HadProjectile::GetDefinition(), G4ParticleTable::GetIon(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4NeutronHPDeExGammas::GetLevel(), G4NeutronHPPhotonDist::GetLevelEnergy(), G4NeutronHPLevel::GetLevelEnergy(), G4NeutronHPDeExGammas::GetLevelEnergy(), G4ReactionProduct::GetMass(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4NucleiProperties::GetNuclearMass(), G4NeutronHPDeExGammas::GetNumberOfLevels(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4NeutronHPPhotonDist::GetPhotons(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), InitDistributionInitialState(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), QI, G4NeutronHPEnAngCorrelation::Sample(), G4NeutronHPEnergyDistribution::Sample(), G4NeutronHPAngular::SampleAndUpdate(), SelectExitChannel(), G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetKineticEnergy(), G4DynamicParticle::SetMomentum(), G4ReactionProduct::SetMomentum(), G4HadFinalState::SetStatusChange(), stopAndKill, theAngularDistribution, G4NeutronHPFinalState::theBaseA, G4NeutronHPFinalState::theBaseZ, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, theGammas, G4NeutronHPFinalState::theResult, and theXsection.

Referenced by G4NeutronHPTInelasticFS::ApplyYourself(), G4NeutronHPPInelasticFS::ApplyYourself(), G4NeutronHPNInelasticFS::ApplyYourself(), G4NeutronHPHe3InelasticFS::ApplyYourself(), G4NeutronHPDInelasticFS::ApplyYourself(), and G4NeutronHPAInelasticFS::ApplyYourself().

00217 {
00218 
00219 // prepare neutron
00220     theResult.Clear();
00221     G4double eKinetic = theTrack.GetKineticEnergy();
00222     const G4HadProjectile *incidentParticle = &theTrack;
00223     G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00224     theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00225     theNeutron.SetKineticEnergy( eKinetic );
00226 
00227 // prepare target
00228     G4int i;
00229     for(i=0; i<50; i++)
00230     { if(theXsection[i] != 0) { break; } } 
00231 
00232     G4double targetMass=0;
00233     G4double eps = 0.0001;
00234     targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
00235                    G4Neutron::Neutron()->GetPDGMass();
00236 //    if(theEnergyAngData[i]!=0)
00237 //        targetMass = theEnergyAngData[i]->GetTargetMass();
00238 //    else if(theAngularDistribution[i]!=0)
00239 //        targetMass = theAngularDistribution[i]->GetTargetMass();
00240 //    else if(theFinalStatePhotons[50]!=0)
00241 //        targetMass = theFinalStatePhotons[50]->GetTargetMass();
00242     G4Nucleus aNucleus;
00243     G4ReactionProduct theTarget; 
00244     G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00245     theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00246 
00247 // prepare the residual mass
00248     G4double residualMass=0;
00249     G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
00250     G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
00251     residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
00252                      G4Neutron::Neutron()->GetPDGMass();
00253 
00254 // prepare energy in target rest frame
00255     G4ReactionProduct boosted;
00256     boosted.Lorentz(theNeutron, theTarget);
00257     eKinetic = boosted.GetKineticEnergy();
00258 //    G4double momentumInCMS = boosted.GetTotalMomentum();
00259   
00260 // select exit channel for composite FS class.
00261     G4int it = SelectExitChannel( eKinetic );
00262    
00263 // set target and neutron in the relevant exit channel
00264     InitDistributionInitialState(theNeutron, theTarget, it);    
00265 
00266     G4ReactionProductVector * thePhotons = 0;
00267     G4ReactionProductVector * theParticles = 0;
00268     G4ReactionProduct aHadron;
00269     aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@    
00270     G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
00271                              (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
00272 //080730c
00273     if ( availableEnergy < 0 )
00274     {
00275        //G4cout << "080730c Adjust availavleEnergy " << G4endl; 
00276        availableEnergy = 0; 
00277     }
00278     G4int nothingWasKnownOnHadron = 0;
00279     G4int dummy;
00280     G4double eGamm = 0;
00281     G4int iLevel=it-1;
00282 
00283 //  TK without photon has it = 0
00284     if( 50 == it ) 
00285     {
00286 
00287 //    TK Excitation level is not determined
00288       iLevel=-1;
00289       aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
00290                                (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
00291 
00292       //aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
00293       //                  std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
00294       //                            aHadron.GetMass()*aHadron.GetMass()));
00295 
00296       //TK add safty 100909
00297       G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
00298       G4double p = 0.0;
00299       if ( p2 > 0.0 ) p = std::sqrt( p ); 
00300 
00301       aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
00302 
00303     }
00304     else
00305     {
00306       while( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; }
00307     }
00308 
00309 
00310     if ( theAngularDistribution[it] != 0 ) // MF4
00311     {
00312       if(theEnergyDistribution[it]!=0) // MF5
00313       {
00314         //************************************************************
00315         /*
00316         aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
00317         G4double eSecN = aHadron.GetKineticEnergy();
00318         */
00319         //************************************************************
00320         //EMendoza --> maximum allowable energy should be taken into account.
00321         G4double dqi = 0.0;
00322         if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
00323         G4double MaxEne=eKinetic+dqi;
00324         G4double eSecN;
00325         do{
00326           eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
00327         }while(eSecN>MaxEne);
00328         aHadron.SetKineticEnergy(eSecN);
00329         //************************************************************
00330         eGamm = eKinetic-eSecN;
00331         for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
00332         {
00333           if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
00334         }
00335         G4double random = 2*G4UniformRand();
00336         iLevel+=G4int(random);
00337         if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
00338       }
00339       else
00340       {
00341         G4double eExcitation = 0;
00342         if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
00343         while (eKinetic-eExcitation < 0 && iLevel>0)
00344         {
00345           iLevel--;
00346           eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
00347         }
00348         //110610TK BEGIN
00349         //Use QI value for calculating excitation energy of residual.
00350         G4bool useQI=false;
00351         G4double dqi = QI[it]; 
00352         if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range
00353  
00354         if ( useQI ) 
00355         {
00356            // QI introudced since G4NDL3.15
00357            eExcitation = -QI[it];
00358            //Re-evluate iLevel based on this eExcitation 
00359            iLevel = 0;
00360            G4bool find = false;
00361            G4int imaxEx = 0;
00362            while( !theGammas.GetLevel(iLevel+1) == 0 ) 
00363            { 
00364               G4double maxEx = 0.0;
00365               if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 
00366               {
00367                  maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();  
00368                  imaxEx = iLevel;
00369               }
00370               if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 
00371               {
00372                  find = true; 
00373                  iLevel--; 
00374                  // very small eExcitation, iLevel becomes -1, this is protected below.
00375                  if ( iLevel == -1 ) iLevel = 0; // But cause energy trouble. 
00376                  break;
00377               }
00378               iLevel++; 
00379            }
00380            // In case, cannot find proper level, then use the maximum level. 
00381            if ( !find ) iLevel = imaxEx;
00382         }
00383         //110610TK END
00384         
00385         if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0) 
00386         {
00387           throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
00388         }
00389         if(eKinetic-eExcitation < 0) eExcitation = 0;
00390         if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
00391         
00392       }
00393       theAngularDistribution[it]->SampleAndUpdate(aHadron);
00394 
00395       if( theFinalStatePhotons[it] == 0 )
00396       {
00397         //G4cout << "110610 USE Gamma Level" << G4endl;
00398 // TK comment Most n,n* eneter to this  
00399         thePhotons = theGammas.GetDecayGammas(iLevel);
00400         eGamm -= theGammas.GetLevelEnergy(iLevel);
00401         if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
00402         {
00403           G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
00404           theRestEnergy->SetDefinition(G4Gamma::Gamma());
00405           theRestEnergy->SetKineticEnergy(eGamm);
00406           G4double costh = 2.*G4UniformRand()-1.;
00407           G4double phi = twopi*G4UniformRand();
00408           theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi), 
00409                                      eGamm*std::sin(std::acos(costh))*std::sin(phi),
00410                                      eGamm*costh);
00411           if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
00412           thePhotons->push_back(theRestEnergy);
00413         }
00414       }
00415     }
00416     else if(theEnergyAngData[it] != 0) // MF6  
00417     {
00418       theParticles = theEnergyAngData[it]->Sample(eKinetic);
00419     }
00420     else
00421     {
00422       // @@@ what to do, if we have photon data, but no info on the hadron itself
00423       nothingWasKnownOnHadron = 1;
00424     }
00425 
00426     //G4cout << "theFinalStatePhotons it " << it << G4endl;
00427     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
00428     //G4cout << "theFinalStatePhotons it " << it << G4endl;
00429     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
00430     //G4cout << "thePhotons " << thePhotons << G4endl;
00431 
00432     if ( theFinalStatePhotons[it] != 0 ) 
00433     {
00434        // the photon distributions are in the Nucleus rest frame.
00435        // TK residual rest frame
00436       G4ReactionProduct boosted_tmp;
00437       boosted_tmp.Lorentz(theNeutron, theTarget);
00438       G4double anEnergy = boosted_tmp.GetKineticEnergy();
00439       thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
00440       G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
00441       G4double testEnergy = 0;
00442       if(thePhotons!=0 && thePhotons->size()!=0)
00443       { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
00444       if(theFinalStatePhotons[it]->NeedsCascade())
00445       {
00446         while(aBaseEnergy>0.01*keV)
00447         {
00448           // cascade down the levels
00449           G4bool foundMatchingLevel = false;
00450           G4int closest = 2;
00451           G4double deltaEold = -1;
00452           for(G4int j=1; j<it; j++)
00453           {
00454             if(theFinalStatePhotons[j]!=0) 
00455             {
00456               testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
00457             }
00458             else
00459             {
00460               testEnergy = 0;
00461             }
00462             G4double deltaE = std::abs(testEnergy-aBaseEnergy);
00463             if(deltaE<0.1*keV)
00464             {
00465               G4ReactionProductVector * theNext = 
00466                 theFinalStatePhotons[j]->GetPhotons(anEnergy);
00467               thePhotons->push_back(theNext->operator[](0));
00468               aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
00469               delete theNext;
00470               foundMatchingLevel = true;
00471               break; // ===>
00472             }
00473             if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
00474             {
00475               closest = j;
00476               deltaEold = deltaE;     
00477             }
00478           } // <=== the break goes here.
00479           if(!foundMatchingLevel)
00480           {
00481             G4ReactionProductVector * theNext = 
00482                theFinalStatePhotons[closest]->GetPhotons(anEnergy);
00483             thePhotons->push_back(theNext->operator[](0));
00484             aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
00485             delete theNext;
00486           }
00487         } 
00488       }
00489     }
00490     unsigned int i0;
00491     if(thePhotons!=0)
00492     {
00493       for(i0=0; i0<thePhotons->size(); i0++)
00494       {
00495         // back to lab
00496         thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
00497       }
00498     }
00499     //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
00500     if(nothingWasKnownOnHadron)
00501     {
00502 //    TKDB 100405
00503 //    In this case, hadron should be isotropic in CM
00504 //    mu and p should be correlated
00505 //
00506       G4double totalPhotonEnergy = 0.0;
00507       if ( thePhotons != 0 )
00508       {
00509          unsigned int nPhotons = thePhotons->size();
00510          unsigned int ii0;
00511          for ( ii0=0; ii0<nPhotons; ii0++)
00512          {
00513             //thePhotons has energies at LAB system 
00514             totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
00515          }
00516       }
00517 
00518       //isotropic distribution in CM 
00519       G4double mu = 1.0 - 2 * G4UniformRand();
00520 
00521       // need momentums in target rest frame;
00522       G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
00523       G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
00524       G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum();
00525 
00526       G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) ); 
00527       G4DynamicParticle* targ = new G4DynamicParticle( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy )  , G4ThreeVector(0) );
00528       G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) );  // will be fill momentum
00529 
00530       two_body_reaction ( proj , targ , hadron , mu );
00531 
00532       G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
00533       G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
00534       aHadron.SetMomentum( hadron_in_LAB.v() );
00535       aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
00536 
00537       delete proj;
00538       delete targ; 
00539       delete hadron;
00540 
00541 //TKDB 100405
00542 /*
00543       G4double totalPhotonEnergy = 0;
00544       if(thePhotons!=0)
00545       {
00546         unsigned int nPhotons = thePhotons->size();
00547         unsigned int i0;
00548         for(i0=0; i0<nPhotons; i0++)
00549         {
00550           totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy();
00551         }
00552       }
00553       availableEnergy -= totalPhotonEnergy;
00554       residualMass += totalPhotonEnergy/G4Neutron::Neutron()->GetPDGMass();
00555       aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
00556                                (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
00557       G4double CosTheta = 1.0 - 2.0*G4UniformRand();
00558       G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
00559       G4double Phi = twopi*G4UniformRand();
00560       G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta);
00561       //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
00562       //                                 aHadron.GetMass()*aHadron.GetMass()));
00563       G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass();
00564 
00565       G4double p = 0.0;
00566       if ( p2 > 0.0 )
00567          p = std::sqrt ( p2 ); 
00568 
00569       aHadron.SetMomentum( Vector*p ); 
00570 */
00571 
00572     }
00573 
00574 // fill the result
00575 // Beware - the recoil is not necessarily in the particles...
00576 // Can be calculated from momentum conservation?
00577 // The idea is that the particles ar emitted forst, and the gammas only once the
00578 // recoil is on the residual; assumption is that gammas do not contribute to 
00579 // the recoil.
00580 // This needs more design @@@
00581 
00582     G4int nSecondaries = 2; // the hadron and the recoil
00583     G4bool needsSeparateRecoil = false;
00584     G4int totalBaryonNumber = 0;
00585     G4int totalCharge = 0;
00586     G4ThreeVector totalMomentum(0);
00587     if(theParticles != 0) 
00588     {
00589       nSecondaries = theParticles->size();
00590       G4ParticleDefinition * aDef;
00591       unsigned int ii0;
00592       for(ii0=0; ii0<theParticles->size(); ii0++)
00593       {
00594         aDef = theParticles->operator[](ii0)->GetDefinition();
00595         totalBaryonNumber+=aDef->GetBaryonNumber();
00596         totalCharge+=G4int(aDef->GetPDGCharge()+eps);
00597         totalMomentum += theParticles->operator[](ii0)->GetMomentum();
00598       } 
00599       if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber())) 
00600       {
00601         needsSeparateRecoil = true;
00602         nSecondaries++;
00603         residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()
00604                           -totalBaryonNumber);
00605         residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge()
00606                           -totalCharge);
00607       }
00608     }
00609     
00610     G4int nPhotons = 0;
00611     if(thePhotons!=0) { nPhotons = thePhotons->size(); }
00612     nSecondaries += nPhotons;
00613         
00614     G4DynamicParticle * theSec;
00615     
00616     if( theParticles==0 )
00617     {
00618       theSec = new G4DynamicParticle;   
00619       theSec->SetDefinition(aHadron.GetDefinition());
00620       theSec->SetMomentum(aHadron.GetMomentum());
00621       theResult.AddSecondary(theSec);    
00622  
00623         aHadron.Lorentz(aHadron, theTarget);
00624         G4ReactionProduct theResidual;   
00625         theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
00626                                   ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
00627         theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
00628 
00629         //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
00630         //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
00631         G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
00632         theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
00633 
00634         theResidual.Lorentz(theResidual, -1.*theTarget);
00635         G4ThreeVector totalPhotonMomentum(0,0,0);
00636         if(thePhotons!=0)
00637         {
00638           for(i=0; i<nPhotons; i++)
00639           {
00640             totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
00641           }
00642         }
00643         theSec = new G4DynamicParticle;   
00644         theSec->SetDefinition(theResidual.GetDefinition());
00645         theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
00646         theResult.AddSecondary(theSec);    
00647     }
00648     else
00649     {
00650       for(i0=0; i0<theParticles->size(); i0++)
00651       {
00652         theSec = new G4DynamicParticle; 
00653         theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
00654         theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
00655         theResult.AddSecondary(theSec); 
00656         delete theParticles->operator[](i0); 
00657       } 
00658       delete theParticles;
00659       if(needsSeparateRecoil && residualZ!=0)
00660       {
00661         G4ReactionProduct theResidual;   
00662         theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
00663                                   ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
00664         G4double resiualKineticEnergy  = theResidual.GetMass()*theResidual.GetMass();
00665                  resiualKineticEnergy += totalMomentum*totalMomentum;
00666                  resiualKineticEnergy  = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
00667 //        cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
00668         theResidual.SetKineticEnergy(resiualKineticEnergy);
00669 
00670         //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
00671         //theResidual.SetMomentum(-1.*totalMomentum);
00672         //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
00673         //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
00674 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
00675         theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
00676 
00677         theSec = new G4DynamicParticle;   
00678         theSec->SetDefinition(theResidual.GetDefinition());
00679         theSec->SetMomentum(theResidual.GetMomentum());
00680         theResult.AddSecondary(theSec);  
00681       }  
00682     }
00683     if(thePhotons!=0)
00684     {
00685       for(i=0; i<nPhotons; i++)
00686       {
00687         theSec = new G4DynamicParticle;    
00688         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00689         //theSec->SetDefinition(G4Gamma::Gamma());
00690         theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
00691         //But never cause real effect at least with G4NDL3.13 TK
00692         theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
00693         theResult.AddSecondary(theSec); 
00694         delete thePhotons->operator[](i);
00695       }
00696 // some garbage collection
00697       delete thePhotons;
00698     }
00699 
00700 //080721 
00701    G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
00702    G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
00703    G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
00704    G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
00705    adjust_final_state ( init_4p_lab ); 
00706 
00707 // clean up the primary neutron
00708     theResult.SetStatusChange(stopAndKill);
00709 }

virtual G4NeutronHPVector* G4NeutronHPInelasticCompFS::GetXsec (  )  [inline, virtual]

Reimplemented from G4NeutronHPFinalState.

Definition at line 83 of file G4NeutronHPInelasticCompFS.hh.

References theXsection.

Referenced by SelectExitChannel().

00083 { return theXsection[50]; }

virtual G4double G4NeutronHPInelasticCompFS::GetXsec ( G4double  anEnergy  )  [inline, virtual]

Reimplemented from G4NeutronHPFinalState.

Definition at line 79 of file G4NeutronHPInelasticCompFS.hh.

References theXsection.

00080   {
00081     return std::max(0., theXsection[50]->GetY(anEnergy));
00082   }

void G4NeutronHPInelasticCompFS::Init ( G4double  A,
G4double  Z,
G4int  M,
G4String dirName,
G4String aSFType 
) [virtual]

Implements G4NeutronHPFinalState.

Reimplemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

Definition at line 78 of file G4NeutronHPInelasticCompFS.cc.

References G4cout, G4endl, gammaPath, G4NeutronHPDataUsed::GetName(), G4NeutronHPNames::GetName(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPEnAngCorrelation::Init(), G4NeutronHPEnergyDistribution::Init(), G4NeutronHPAngular::Init(), G4NeutronHPVector::Init(), G4NeutronHPPhotonDist::InitAngular(), G4NeutronHPPhotonDist::InitEnergies(), G4NeutronHPPhotonDist::InitMean(), G4NeutronHPPhotonDist::InitPartials(), LR, QI, G4NeutronHPFinalState::SetAZMs(), theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, G4NeutronHPFinalState::theNames, G4NeutronHPFinalState::theNDLDataA, G4NeutronHPFinalState::theNDLDataZ, and theXsection.

Referenced by G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), G4NeutronHPDInelasticFS::Init(), and G4NeutronHPAInelasticFS::Init().

00079 {
00080   gammaPath = "/Inelastic/Gammas/";
00081     if(!getenv("G4NEUTRONHPDATA")) 
00082        throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00083   G4String tBase = getenv("G4NEUTRONHPDATA");
00084   gammaPath = tBase+gammaPath;
00085   G4String tString = dirName;
00086   G4bool dbool;
00087   G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
00088   G4String filename = aFile.GetName();
00089   SetAZMs( A, Z, M, aFile ); 
00090   //theBaseA = aFile.GetA();
00091   //theBaseZ = aFile.GetZ();
00092   //theNDLDataA = (int)aFile.GetA();
00093   //theNDLDataZ = aFile.GetZ();
00094   //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
00095   if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
00096   {
00097     if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
00098     hasAnyData = false;
00099     hasFSData = false; 
00100     hasXsec = false;
00101     return;
00102   }
00103   //theBaseA = A;
00104   //theBaseZ = G4int(Z+.5);
00105   std::ifstream theData(filename, std::ios::in);
00106   if(!theData)
00107   {
00108     hasAnyData = false;
00109     hasFSData = false; 
00110     hasXsec = false;
00111     theData.close();
00112     return;
00113   }
00114   // here we go
00115   G4int infoType, dataType, dummy;
00116   G4int sfType, it;
00117   hasFSData = false; 
00118   while (theData >> infoType)
00119   {
00120     hasFSData = true; 
00121     theData >> dataType;
00122     theData >> sfType >> dummy;
00123     it = 50;
00124     if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
00125     if(dataType==3) 
00126     {
00127       //theData >> dummy >> dummy;
00128       //TK110430
00129       // QI and LR introudced since G4NDL3.15
00130       G4double dqi;
00131       G4int ilr;
00132       theData >> dqi >> ilr;
00133 
00134       QI[ it ] = dqi*eV;
00135       LR[ it ] = ilr;
00136       theXsection[it] = new G4NeutronHPVector;
00137       G4int total;
00138       theData >> total;
00139       theXsection[it]->Init(theData, total, eV);
00140       //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
00141     }
00142     else if(dataType==4)
00143     {
00144       theAngularDistribution[it] = new G4NeutronHPAngular;
00145       theAngularDistribution[it]->Init(theData);
00146     }
00147     else if(dataType==5)
00148     {
00149       theEnergyDistribution[it] = new G4NeutronHPEnergyDistribution;
00150       theEnergyDistribution[it]->Init(theData); 
00151     }
00152     else if(dataType==6)
00153     {
00154       theEnergyAngData[it] = new G4NeutronHPEnAngCorrelation;
00155       theEnergyAngData[it]->Init(theData);
00156     }
00157     else if(dataType==12)
00158     {
00159       theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
00160       theFinalStatePhotons[it]->InitMean(theData);
00161     }
00162     else if(dataType==13)
00163     {
00164       theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
00165       theFinalStatePhotons[it]->InitPartials(theData);
00166     }
00167     else if(dataType==14)
00168     {
00169       theFinalStatePhotons[it]->InitAngular(theData);
00170     }
00171     else if(dataType==15)
00172     {
00173       theFinalStatePhotons[it]->InitEnergies(theData);
00174     }
00175     else
00176     {
00177       throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS");
00178     }
00179   }
00180   theData.close();
00181 }

void G4NeutronHPInelasticCompFS::InitDistributionInitialState ( G4ReactionProduct aNeutron,
G4ReactionProduct aTarget,
G4int  it 
) [inline]

Definition at line 86 of file G4NeutronHPInelasticCompFS.hh.

References G4NeutronHPEnAngCorrelation::SetNeutron(), G4NeutronHPAngular::SetNeutron(), G4NeutronHPEnAngCorrelation::SetTarget(), G4NeutronHPAngular::SetTarget(), theAngularDistribution, and theEnergyAngData.

Referenced by CompositeApply().

00089   {
00090     if(theAngularDistribution[it]!=0) 
00091     {
00092       theAngularDistribution[it]->SetTarget(aTarget);
00093       theAngularDistribution[it]->SetNeutron(aNeutron);
00094     }
00095     if(theEnergyAngData[it]!=0)
00096     {
00097       theEnergyAngData[it]->SetTarget(aTarget);
00098       theEnergyAngData[it]->SetNeutron(aNeutron);
00099     }
00100   }

void G4NeutronHPInelasticCompFS::InitGammas ( G4double  AR,
G4double  ZR 
)

Definition at line 57 of file G4NeutronHPInelasticCompFS.cc.

References gammaPath, G4NeutronHPDeExGammas::Init(), and theGammas.

Referenced by G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), G4NeutronHPDInelasticFS::Init(), and G4NeutronHPAInelasticFS::Init().

00058 {
00059   //   char the[100] = {""};
00060   //   std::ostrstream ost(the, 100, std::ios::out);
00061   //   ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
00062   //   G4String * aName = new G4String(the);
00063   //   std::ifstream from(*aName, std::ios::in);
00064 
00065    std::ostringstream ost;
00066    ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
00067    G4String aName = ost.str();
00068    std::ifstream from(aName, std::ios::in);
00069 
00070    if(!from) return; // no data found for this isotope
00071    //   std::ifstream theGammaData(*aName, std::ios::in);
00072    std::ifstream theGammaData(aName, std::ios::in);
00073     
00074    theGammas.Init(theGammaData);
00075    //   delete aName;
00076 }

virtual G4NeutronHPFinalState* G4NeutronHPInelasticCompFS::New (  )  [pure virtual]

Implements G4NeutronHPFinalState.

Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

G4int G4NeutronHPInelasticCompFS::SelectExitChannel ( G4double  eKinetic  ) 

Definition at line 183 of file G4NeutronHPInelasticCompFS.cc.

References G4UniformRand, GetXsec(), and theXsection.

Referenced by CompositeApply().

00184 {
00185 
00186 // it = 0 has without Photon
00187   G4double running[50];
00188   running[0] = 0;
00189   unsigned int i;
00190   for(i=0; i<50; i++)
00191   {
00192     if(i!=0) running[i]=running[i-1];
00193     if(theXsection[i] != 0) 
00194     {
00195       running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
00196     }
00197   }
00198   G4double random = G4UniformRand();
00199   G4double sum = running[49];
00200   G4int it = 50;
00201   if(0!=sum)
00202   {
00203     G4int i0;
00204     for(i0=0; i0<50; i0++)
00205     {
00206       it = i0;
00207       if(random < running[i0]/sum) break;
00208     }
00209   }
00210 //debug:  it = 1;
00211   return it;
00212 }


Field Documentation

G4String G4NeutronHPInelasticCompFS::gammaPath [protected]

Definition at line 112 of file G4NeutronHPInelasticCompFS.hh.

Referenced by Init(), and InitGammas().

std::vector<G4int > G4NeutronHPInelasticCompFS::LR [protected]

Definition at line 119 of file G4NeutronHPInelasticCompFS.hh.

Referenced by G4NeutronHPInelasticCompFS(), and Init().

std::vector< G4double > G4NeutronHPInelasticCompFS::QI [protected]

Definition at line 118 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), and Init().

G4NeutronHPAngular* G4NeutronHPInelasticCompFS::theAngularDistribution[51] [protected]

Definition at line 106 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4NeutronHPInelasticCompFS().

G4double G4NeutronHPInelasticCompFS::theCurrentA [protected]

Definition at line 114 of file G4NeutronHPInelasticCompFS.hh.

G4double G4NeutronHPInelasticCompFS::theCurrentZ [protected]

Definition at line 115 of file G4NeutronHPInelasticCompFS.hh.

G4NeutronHPEnAngCorrelation* G4NeutronHPInelasticCompFS::theEnergyAngData[51] [protected]

Definition at line 107 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4NeutronHPInelasticCompFS().

G4NeutronHPEnergyDistribution* G4NeutronHPInelasticCompFS::theEnergyDistribution[51] [protected]

Definition at line 105 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), and ~G4NeutronHPInelasticCompFS().

G4NeutronHPPhotonDist* G4NeutronHPInelasticCompFS::theFinalStatePhotons[51] [protected]

Definition at line 109 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), and ~G4NeutronHPInelasticCompFS().

G4NeutronHPDeExGammas G4NeutronHPInelasticCompFS::theGammas [protected]

Definition at line 111 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), and InitGammas().

G4NeutronHPVector* G4NeutronHPInelasticCompFS::theXsection[51] [protected]

Definition at line 104 of file G4NeutronHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), GetXsec(), Init(), SelectExitChannel(), and ~G4NeutronHPInelasticCompFS().


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