G4NeutronHPInelasticCompFS.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 // neutron_hp -- source file
00027 // J.P. Wellisch, Nov-1996
00028 // A prototype of the low energy neutron transport model.
00029 //
00030 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
00031 // 070606 bug fix and migrate to enable to Partial cases by T. Koi 
00032 // 080603 bug fix for Hadron Hyper News #932 by T. Koi 
00033 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
00034 // 080717 bug fix of calculation of residual momentum by T. Koi
00035 // 080801 protect negative avalable energy by T. Koi
00036 //        introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
00037 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
00038 // 090514 Fix bug in IC electron emission case 
00039 //        Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
00040 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM 
00041 //        add two_body_reaction
00042 // 100909 add safty 
00043 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation  
00044 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
00045 //
00046 #include "G4NeutronHPInelasticCompFS.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4Nucleus.hh"
00050 #include "G4NucleiProperties.hh"
00051 #include "G4He3.hh"
00052 #include "G4Alpha.hh"
00053 #include "G4Electron.hh"
00054 #include "G4NeutronHPDataUsed.hh"
00055 #include "G4ParticleTable.hh"
00056 
00057 void G4NeutronHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
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 }
00077 
00078 void G4NeutronHPInelasticCompFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType)
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 }
00182 
00183 G4int G4NeutronHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
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 }
00213 
00214 
00215                                                                                                        //n,p,d,t,he3,a
00216 void G4NeutronHPInelasticCompFS::CompositeApply(const G4HadProjectile & theTrack, G4ParticleDefinition * aDefinition)
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 }
00710 
00711 
00712 
00713 #include "G4RotationMatrix.hh" 
00714 void G4NeutronHPInelasticCompFS::two_body_reaction ( G4DynamicParticle* proj, G4DynamicParticle* targ, G4DynamicParticle* hadron, G4double mu ) 
00715 {
00716 
00717 // Target rest flame
00718 // 4vector in targ rest frame;
00719 // targ could have excitation energy (photon energy will be emiited) tricky but,,,
00720 
00721    G4LorentzVector before = proj->Get4Momentum() + targ->Get4Momentum();
00722 
00723    G4ThreeVector p3_proj = proj->GetMomentum();
00724    G4ThreeVector d = p3_proj.unit();
00725    G4RotationMatrix rot; 
00726    G4RotationMatrix rot1; 
00727    rot1.setPhi( pi/2 + d.phi() );
00728    G4RotationMatrix rot2; 
00729    rot2.setTheta( d.theta() );
00730    rot=rot2*rot1;
00731    proj->SetMomentum( rot*p3_proj );
00732 
00733 // Now proj only has pz component;
00734 
00735 // mu in CM system 
00736 
00737    //Valid only for neutron incidence
00738    G4DynamicParticle* residual = new G4DynamicParticle ( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)( targ->GetDefinition()->GetPDGCharge() - hadron->GetDefinition()->GetPDGCharge() ) , (G4int)(targ->GetDefinition()->GetBaryonNumber() - hadron->GetDefinition()->GetBaryonNumber()+1) , 0 ) , G4ThreeVector(0) ); 
00739 
00740    G4double Q = proj->GetDefinition()->GetPDGMass() + targ->GetDefinition()->GetPDGMass() 
00741               - ( hadron->GetDefinition()->GetPDGMass() + residual->GetDefinition()->GetPDGMass() );
00742 
00743    // Non Relativistic Case 
00744    G4double A = targ->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
00745    G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass(); 
00746    G4double E1 = proj->GetKineticEnergy();
00747 
00748 // 101111
00749 // In _nat_ data (Q+E1) could become negative value, following line is safty for this case.
00750    //if ( (Q+E1) < 0 ) 
00751    if ( ( 1 + (1+A)/A*Q/E1 ) < 0 ) 
00752    {
00753 // 1.0e-6 eV is additional safty for numeric precision 
00754       Q = -( A/(1+A)*E1 ) + 1.0e-6*eV;
00755    }
00756 
00757    G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
00758    G4double gamma = AA/(A+1-AA)*beta;
00759    G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
00760    G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
00761    if ( omega3 > 1.0 ) omega3 = 1.0;
00762 
00763    G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
00764    G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
00765    if ( omega4 > 1.0 ) omega4 = 1.0;
00766 
00767    hadron->SetKineticEnergy ( E3 );
00768    
00769    G4double M = hadron->GetDefinition()->GetPDGMass();
00770    G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
00771    G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
00772 
00773    G4double M4 = residual->GetDefinition()->GetPDGMass();
00774    G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
00775    G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
00776 
00777 // Rotate to orginal target rest flame.
00778    p *= rot.inverse();
00779    hadron->SetMomentum( p );
00780 // Now hadron had 4 momentum in target rest flame 
00781 
00782 // TypeA
00783    p4 *= rot.inverse();
00784    residual->SetMomentum ( p4 );
00785 
00786 //TypeB1
00787    //residual->Set4Momentum ( p4_residual );
00788 //TypeB2
00789    //residual->SetMomentum ( p4_residual.v() );
00790 
00791 // Type A make difference in Momenutum
00792 // Type B1 make difference in Mass of residual
00793 // Type B2 make difference in total energy.
00794 
00795    delete residual;
00796 
00797 }

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