| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

DsG4NeutronHPCapture Class Reference

#include <DsG4NeutronHPCapture.h>

Collaboration diagram for DsG4NeutronHPCapture:

[legend]
List of all members.

Public Member Functions

 DsG4NeutronHPCapture ()
virtual ~DsG4NeutronHPCapture ()
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
void passNeutronCaptureInfoTool (INeutronCaptureInfo *)

Public Attributes

G4bool DebugMe

Private Attributes

G4double * xSec
G4NeutronHPChannel * theCapture
G4String dirName
G4int numEle
G4int it
G4HadFinalState theResult
G4HadFinalState * result
INeutronCaptureInfo * m_capinfo
G4DhNeutronCapture nCapture

Detailed Description

Definition at line 24 of file DsG4NeutronHPCapture.h.

Constructor & Destructor Documentation

DsG4NeutronHPCapture::DsG4NeutronHPCapture (  ) 

Definition at line 32 of file DsG4NeutronHPCapture.cc.

00033 {
00034   DebugMe = false;
00035   if (DebugMe) G4cout  << "Begin DsG4NeutronCapture constructor " << G4endl;
00036 
00037   SetMinEnergy( 0.0 );
00038   SetMaxEnergy( 20.*MeV );
00039   
00040   char* envstr = getenv("G4NEUTRONHPDATA");
00041 
00042   if(!envstr) {
00043     throw G4HadronicException(__FILE__, __LINE__, 
00044                               "Please setenv G4NEUTRONHPDATA "
00045                               "to point to the neutron cross-section files.");
00046   }
00047   G4String tString = "/Capture/";
00048   dirName = envstr;
00049   dirName = dirName + tString;
00050   numEle = G4Element::GetNumberOfElements();
00051 
00052   G4cout << "Initializing DsG4NeutronCapture for " << numEle << " neutron HP channels" << G4endl;
00053   
00054   theCapture = new G4NeutronHPChannel[numEle];
00055 
00056   G4NeutronHPCaptureFS * theFS = new G4NeutronHPCaptureFS; 
00057   for (G4int i=0; i<numEle; i++) {
00058     
00059     if (DebugMe) G4cout << "initializing theCapture "<<i<<" "<< numEle
00060             << " name " << (*(G4Element::GetElementTable()))[i]->GetName()
00061             << G4endl; 
00062     // DEBUG:
00063     
00064     G4cout << "initializing theCapture "<<i<<" "<< numEle
00065          << ", name " << (*(G4Element::GetElementTable()))[i]->GetName() 
00066          << ", symbol " << (*(G4Element::GetElementTable()))[i]->GetSymbol() 
00067          << ", z   " << (*(G4Element::GetElementTable()))[i]->GetZ() 
00068          << ", n   " << (*(G4Element::GetElementTable()))[i]->GetN() 
00069          << ", a   " << (*(G4Element::GetElementTable()))[i]->GetA() 
00070          << ", nat " << (*(G4Element::GetElementTable()))[i]->GetNaturalAbandancesFlag() 
00071          << ", Niso   " << (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() 
00072          << G4endl;
00073 
00074     // Cast Z to nearest integer
00075     int elemZ = (int)((*(G4Element::GetElementTable()))[i]->GetZ() + 0.5);
00076 
00077     // To use new NNDC Gd Capture Gammas, change this to true
00078     bool useNNDC_GdCapture = false;
00079 
00080     if(elemZ == 64 && !useNNDC_GdCapture) { // Gd.
00081       // for Gd, DsG4GdNeutronHPCaptureFS is invoked.
00082       DsG4GdNeutronHPCaptureFS * theGdFS = new DsG4GdNeutronHPCaptureFS;
00083       theCapture[i].Init((*(G4Element::GetElementTable()))[i], dirName);
00084       if (DebugMe) G4cout << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
00085       theCapture[i].Register(theGdFS);
00086       delete theGdFS;
00087     }
00088     else if(elemZ == 6 /*C*/
00089             || elemZ == 7 /*N*/
00090             || elemZ == 8 /*O*/
00091             || elemZ == 14 /*Si*/
00092             || elemZ == 15 /*P*/
00093             || elemZ == 16 /*S*/
00094             || elemZ == 24 /*Cr*/
00095             || elemZ == 25 /*Mn*/
00096             || elemZ == 26 /*Fe*/
00097             || elemZ == 28 /*Ni*/
00098             || (elemZ == 64 && useNNDC_GdCapture) /*Gd*/){
00099       // For these elements, we use hand-generated tables of correlated gammas.
00100       // These tables ensure conservation of energy
00101       theCapture[i].Init((*(G4Element::GetElementTable()))[i], dirName);
00102       DsG4NNDCNeutronHPCaptureFS * theNNDCFS = new DsG4NNDCNeutronHPCaptureFS();
00103       theCapture[i].Register(theNNDCFS);
00104       delete theNNDCFS;
00105     }
00106     else { 
00107       theCapture[i].Init((*(G4Element::GetElementTable()))[i], dirName);
00108       theCapture[i].Register(theFS);
00109     }
00110   }
00111   delete theFS;
00112 }

DsG4NeutronHPCapture::~DsG4NeutronHPCapture (  )  [virtual]

Definition at line 114 of file DsG4NeutronHPCapture.cc.

00115 {
00116   delete [] theCapture;
00117 }

Member Function Documentation

G4HadFinalState * DsG4NeutronHPCapture::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  aTargetNucleus 
)

Definition at line 121 of file DsG4NeutronHPCapture.cc.

00123 {
00124   result = new G4HadFinalState(); 
00125 
00126   // Initialize
00127   G4int gammaNum = 0; 
00128   G4double capGammaE[20] = {0}; 
00129   G4double capTime = 0; 
00130   G4double capGammaEsum = 0;
00131   
00132   //if(getenv("NeutronHPCapture")) 
00133   if ( DebugMe ) {G4cout <<" ### DsG4NeutronHPCapture called"<< G4endl;} 
00134 
00135   // get cross-sections for elements in current material
00136   const G4Material * theMaterial = aTrack.GetMaterial();
00137   G4int n = theMaterial->GetNumberOfElements();
00138   G4int index = theMaterial->GetElement(0)->GetIndex();
00139   if(n!=1) {
00140     xSec = new G4double[n];
00141     G4double sum=0;
00142     G4int i;
00143     const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
00144     G4double rWeight;    
00145     G4NeutronHPThermalBoost aThermalE;
00146     for (i=0; i<n; i++) {
00147       index = theMaterial->GetElement(i)->GetIndex();
00148       rWeight = NumAtomsPerVolume[i];
00149       xSec[i] = theCapture[index].GetXsec(
00150                                           aThermalE.GetThermalEnergy(aTrack,
00151                                                                      theMaterial->GetElement(i),
00152                                                                      theMaterial->GetTemperature()));
00153 
00154       xSec[i] *= rWeight;
00155       sum+=xSec[i];
00156     }
00157     
00158     // determine the target nucleus
00159     G4double random = G4UniformRand();
00160     G4double running = 0;
00161     for (i=0; i<n; i++) {
00162       running += xSec[i];
00163       index = theMaterial->GetElement(i)->GetIndex();
00164       if(random<=running/sum) break;
00165     }
00166     if(i==n) i=std::max(0, n-1);
00167     delete [] xSec;
00168   }
00169 
00170   capTime = aTrack.GetGlobalTime()/(1000*ns);
00171   
00172   G4String targetname = (*(G4Element::GetElementTable()))[index]->GetName();
00173   G4double targetZ = (*(G4Element::GetElementTable()))[index]->GetZ();
00174   G4double targetA = (*(G4Element::GetElementTable()))[index]->GetN(); 
00175   
00176   result = new G4HadFinalState();
00177 
00178   // Allow up to 100 tries to get physically meaningful capture gamma energy and number of gammas
00179   G4int tries = 100;  
00180   G4bool done = false;  
00181   while (!done) { 
00182 
00183     gammaNum = 0; 
00184     capGammaEsum = 0;
00185 
00186     result = theCapture[index].ApplyYourself(aTrack);
00187       
00188     G4int num = result->GetNumberOfSecondaries();
00189 
00190     if (DebugMe) G4cout << "DDEBUG: number of secondaries: " << num << G4endl;
00191       
00192     G4String secname;
00193     G4double seckine;
00194     for(int ii=0;ii<num;ii++) {
00195       secname=(result->GetSecondary(ii))->GetParticle()->GetDefinition()->GetParticleName();
00196       seckine=(result->GetSecondary(ii))->GetParticle()->GetKineticEnergy()/MeV;  
00197 
00198       if (DebugMe) G4cout << "   DDEBUG   name: " << secname << G4endl; 
00199       // Attention: the recoiling target is one of the secondaries.
00200       if(secname == "gamma") {
00201         gammaNum++;
00202         capGammaE[gammaNum-1] = seckine;
00203         nCapture.PushCapGammaE(seckine); // Push info into n-cap tool
00204       }
00205       if(secname == "gamma") {
00206         capGammaEsum += seckine;
00207       }    
00208     }//end of loop over secondaries
00209       
00210     // djaffe: decide if gammas are physically reasonable
00211       //         Captures on H and Gd give valid results.
00212       //         Enforce # of gammas and total energy for captures on carbon.
00213       //         For captures on other nuclei, just avoid zero or very small gamma energy.
00214       //         This needs to be fixed in the future.
00215       int iz = (int)(targetZ + 0.5); // convert to integer
00216 
00217       switch (iz) {
00218       case 1: // hydrogen
00219         done = true ; // nH always OK
00220         break;
00221       case 64: // Gd
00222         done = true ; // nGd always OK
00223         break;
00224         //case 6: // carbon
00225         //      done = ((gammaNum == 1) || (gammaNum == 2)) && std::abs(capGammaEsum-4.946*MeV)<0.1*MeV ;
00226         //      break;
00227       default: // everything else
00228         if ( capGammaEsum > 0.01*MeV ) {
00229           G4bool zero = false;
00230           for (int ii=0; ii<gammaNum; ii++) { if (capGammaE[ii] < 0.001*MeV) zero = true; }
00231           done = !zero;
00232         }
00233         break;
00234       }
00235       --tries; // don't try forever
00236       if (tries == 0)   {
00237         done = true; 
00238         G4cout << " DsG4NeutronHPCapture: GIVING UP. Could not achieve acceptable final state gammas  " << G4endl;
00239         G4cout << " DsG4NeutronHPCapture: Z,A " <<  targetZ <<","<<targetA
00240                  <<" "<<targetname << " N(gamma)="
00241                  << gammaNum << " E(gamma)= " ;                             
00242         for (int ii=0;ii<gammaNum;ii++){ G4cout << capGammaE[ii] <<", " ;} 
00243         G4cout << G4endl; 
00244       }
00245     } // !done
00246     
00247     
00248     /* recording the capture information --- Wei Wang, Aug 14, 2008 */
00249     nCapture.SetCapTargetZ(targetZ);
00250     nCapture.SetCapTargetA(targetA);
00251     nCapture.SetCapTime(capTime);
00252     nCapture.SetCapGammaN(gammaNum);
00253     
00254     m_capinfo->addCapture(nCapture);
00255     /* capture information recorded */
00256 
00257     return result; 
00258 }

void DsG4NeutronHPCapture::passNeutronCaptureInfoTool ( INeutronCaptureInfo *   ) 

Definition at line 260 of file DsG4NeutronHPCapture.cc.

00261 {
00262   m_capinfo = p_capinfo;
00263 }

Member Data Documentation

G4bool DsG4NeutronHPCapture::DebugMe

Definition at line 34 of file DsG4NeutronHPCapture.h.

G4double* DsG4NeutronHPCapture::xSec [private]

Definition at line 40 of file DsG4NeutronHPCapture.h.

G4NeutronHPChannel* DsG4NeutronHPCapture::theCapture [private]

Definition at line 41 of file DsG4NeutronHPCapture.h.

G4String DsG4NeutronHPCapture::dirName [private]

Definition at line 42 of file DsG4NeutronHPCapture.h.

G4int DsG4NeutronHPCapture::numEle [private]

Definition at line 43 of file DsG4NeutronHPCapture.h.

G4int DsG4NeutronHPCapture::it [private]

Definition at line 44 of file DsG4NeutronHPCapture.h.

G4HadFinalState DsG4NeutronHPCapture::theResult [private]

Definition at line 46 of file DsG4NeutronHPCapture.h.

G4HadFinalState* DsG4NeutronHPCapture::result [private]

Definition at line 47 of file DsG4NeutronHPCapture.h.

INeutronCaptureInfo* DsG4NeutronHPCapture::m_capinfo [private]

Definition at line 52 of file DsG4NeutronHPCapture.h.

G4DhNeutronCapture DsG4NeutronHPCapture::nCapture [private]

Definition at line 53 of file DsG4NeutronHPCapture.h.

The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |
Generated on Mon Apr 11 20:53:27 2011 for DetSim by doxygen 1.4.7