Message: Re: G4Exception : 001 issued by : G4HadronicProcess Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Idea Re: G4Exception : 001 issued by : G4HadronicProcess 

Keywords: G4Exception : 001 issued by : G4HadronicProcess
Forum: Hadronic Processes
Re: None G4Exception : 001 issued by : G4HadronicProcess (Lam YiHua)
Re: Feedback Re: G4Exception : 001 issued by : G4HadronicProcess (Gunter Folger)
Date: 09 Apr, 2007
From: LAM YiHua <LamYiHua@gmail.com>

Dear Gunter, Dennis and Experts in Hadronic physics,

 

I define elements by means of using isotopes and their natural abundance. But I will face the problems which I had posted at:

 

http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/hadronprocess/672.html

http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/hadronprocess/628.html

 

 

I found that when I modify the original method G4HadronicProcess::ChooseAandZ by discarding the implementation of "static G4StableIsotopes theIso",

which will redefine elements (and may disregard the percentage of components of elements defined by application developers in the instance UserDetectorConstruction),

and replacing them with the code below, both segmentation faults will not appear.

Could you please help me to check whether these codes are legal to use?

 

Many thanks... ;)

 

 

 

 

best regards,

 

YiHua

 

 

_____________________________________________________________

 

G4Element * MyHadronicProcess::ChooseAandZ(const G4DynamicParticle *aParticle, const G4Material *aMaterial )

{ //YiHua altered Jan 2007

  static G4bool noIsotopeWiseCrossSections=getenv("GHAD_DISABLE_ISOTOPE_WISE_CROSS_SECTIONS");

 

  currentZ = 0;

  currentN = 0;

  const G4int                  numberOfElements = aMaterial->GetNumberOfElements();

  const G4ElementVector *theElementVector = aMaterial->GetElementVector();

 

//---------------     ---------------     ---------------     ---------------   

//  First Loop to facilitate single element, single isotope or multiple isotopes

//---------------     ---------------     ---------------     ---------------

if( numberOfElements == 1 )

{

currentZ = G4double( ((*theElementVector)[0])->GetZ());

const G4IsotopeVector *theIsotopeVector = (*theElementVector)[0]->GetIsotopeVector();

 

  if(noIsotopeWiseCrossSections) { currentN = (*theElementVector)[0]->GetN(); }

  else

  {

     G4int NumberOfIsotopes = (*theElementVector)[0]- >GetNumberOfIsotopes();

 

// +++++++++++++++++++++++++++++++++++

// Single Element, Single Isotope

// +++++++++++++++++++++++++++++++++++

  if (NumberOfIsotopes == 1)

  {

     currentN = (*theIsotopeVector)[0]->GetN();

     targetNucleus.SetParameters(currentN, currentZ);

     return (*theElementVector)[0];

  }

// +++++++++++++++++++++++++++++++++++

// Single Element, Multiple Isotopes

// +++++++++++++++++++++++++++++++++++

  else if (NumberOfIsotopes > 1)

  {

     G4double                           *running = new G4double [(*theElementVector)[0]->GetNumberOfIsotopes()];

     G4double *RelativeAbundanceVector = (*theElementVector)[0]- >GetRelativeAbundanceVector();

 

     for (G4int i=0; i<NumberOfIsotopes; i++)

     {

         G4double fracInPercent = RelativeAbundanceVector[i];

         G4double       runningA = (*theIsotopeVector)[i]->GetN(); //get number of nucleon

                             running[i] = fracInPercent*std::pow(runningA, 2./3.);

                 if(i!=0){running[i] += running[i-1];}

     }

 

     G4double trial = G4UniformRand();

     for(G4int k=0; k<NumberOfIsotopes; k++)

     {

         currentN = (*theIsotopeVector)[k]->GetN();

         if(running[k]/running[NumberOfI sotopes-1]>trial) { break;}

     }

     delete [] running;

  }

}

 

// Check for Z > 92)

  if (currentZ > 92) G4Exception("MyHadronicProcess ", "008", FatalException, "Z > 92 not allowed");

 

  targetNucleus.SetParameters(currentN, currentZ);

  return (*theElementVector)[0];

}

 

//---------------     ---------------     ---------------     ---------------

//  Second Loop to facilitate Multi elements, Single isotope or Multi-Isotopes

//---------------     ---------------     ---------------     ---------------

 

const G4double *theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();

G4double                                       aTemp = aMaterial->GetTemperature();

G4double                       crossSectionTotal = 0;

 

std::vector<G4double> runningSum;

 

//

// Calculate crossSectionTotal

//

 

for(G4int i=0; i < numberOfElements; ++i )

{

  runningSum.push_back (theAtomicNumberDensity[i]*dispatch->GetMicroscopicCrossSection( aParticle, (*theElementVector)[i], aTemp));

  crossSectionTotal+=runningSum[i];

}

 

 

//

// Monte Carlo Selection of an element from various elements

//

G4double random = G4UniformRand();

for(G4int i=0; i < numberOfElements; ++i )

{

  if(i!=0){runningSum[i]+=runningSum[i-1];}

 

  if( random<=runningSum[i]/crossSectionTotal )

  {

     currentZ = G4double( ((*theElementVector)[i])->GetZ());

   

     if(noIsotopeWiseCrossSections)

     {  currentN = ((*theElementVector)[i])->GetN();  }

//

// After selection of an element...

// Monte Carlo Selection of an isotope from various isotopes

//

     else

     {

        G4double          *running = new G4double [(*theElementVector)[i]- >GetNumberOfIsotopes()];

        G4int NumberOfIsotopes = (*theElementVector)[i]->GetNumberOfIsotopes();

 

        const G4IsotopeVector  *theIsotopeVector = (*theElementVector)[i]- >GetIsotopeVector();

        G4double         *RelativeAbundanceVector = (*theElementVector)[i]->GetRelativeAbundanceVector();

 

        if (NumberOfIsotopes == 1)

        {

           currentN = G4double( ((*theElementVector)[i])->GetN());

           break;

        }

        else if (NumberOfIsotopes > 1)

        {

           for (G4int j=0; j<NumberOfIsotopes; j++)

           {

              G4double  fracInPercent = RelativeAbundanceVector[j];

              G4double        runningA = (*theIsotopeVector)[j]->GetN(); //get number of nucleon

                                   running[j] = fracInPercent*std::pow(runningA, 2./3.);

                        if(j!=0) running[j] += running[j-1];

           }

         

           G4double trial = G4UniformRand();

           for(G4int k=0; k<NumberOfIsotopes; k++)

           {

              currentN = (*theIsotopeVector)[k]->GetN();

 

              if(running[k]/running[NumberOfIsotopes-1]>trial) { break; }

           }

        }

        delete [] running;

     }

 

     targetNucleus.SetParameters(currentN, currentZ);

     return (*theElementVector)[i];

  }

}

//--------------- --------------- --------------- --------------- ---------------

// Third Loop to facilitate the Last registered element, Single isotope or Multi-Isotopes

//--------------- --------------- --------------- --------------- ---------------

 

currentZ = G4double((*theElementVector)[numberOfElements-1]->GetZ());

 

if(noIsotopeWiseCrossSections)

{

  currentN = (*theElementVector)[numberOfElements-1]->GetN();

}

 

else

{

  G4double                                *running = new G4double [(*theElementVector)[numberOfElements-1]->GetNumberOfIsotopes()];

  G4int                       NumberOfIsotopes = (*theElementVector)[numberOfEle ments-1]->GetNumberOfIsotopes();

const G4IsotopeVector *theIsotopeVector = (*theElementVector)[numberOfElements-1]->GetIsotopeVector();

G4double        *RelativeAbundanceVector = (*theElementVector)[numberOfElements-1]->GetRelativeAbundanceVector();

 

for (G4int i=0; i<NumberOfIsotopes; i++)

{

  G4double fracInPercent = RelativeAbundanceVector[i];

  G4double       runningA = (*theIsotopeVector)[i]->GetN();

                      running[i] = fracInPercent*std::pow(runningA , 2./3.);

 

  if(i!=0){ running[i] += running[i-1];}

}

 

G4double trial = G4UniformRand();

 

for(G4int i=0; i<NumberOfIsotopes; i++)

{

  currentN = (*theIsotopeVector)[i]->GetN();

  if(running[i]/running[NumberOfIsotopes-1]>trial) {break;}

}

delete [] running;

targetNucleus.SetParameters (currentN, currentZ);

return (*theElementVector)[numberOfElements-1];

}

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

 Add Message Add Message
to: "Re: G4Exception : 001 issued by : G4HadronicProcess"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]