Changeset 1055 for trunk/source/processes/electromagnetic/utils/src/G4EmElementSelector.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/electromagnetic/utils/src/G4EmElementSelector.cc (modified) (7 diffs)

trunk/source/processes/electromagnetic/utils/src/G4EmElementSelector.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmElementSelector.cc,v 1.4 2008/08/21 18:53:32 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4EmElementSelector.cc,v 1.10 2009/05/26 16:59:35 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -58 +58 @@
                                          G4double emin, 
                                          G4double emax,
-                                         G4bool spline):
+                                         G4bool /*spline*/):
   model(mod), material(mat), nbins(bins), cutEnergy(-1.0), 
   lowEnergy(emin), highEnergy(emax)
@@ -65 +65 @@
   nElmMinusOne = n - 1;
   theElementVector = material->GetElementVector();
+  element = (*theElementVector)[0];
   if(nElmMinusOne > 0) {
-    for(G4int i=0; i<nElmMinusOne; i++) {
+    xSections.reserve(n);
+    for(G4int i=0; i<n; i++) {
       G4PhysicsLogVector* v = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins);
-      v->SetSpline(spline);
+      //v->SetSpline(spline);
       xSections.push_back(v);
     }
@@ -80 +82 @@
 {
   if(nElmMinusOne > 0) {
-    for(G4int i=0; i<nElmMinusOne; i++) {
+    for(G4int i=0; i<=nElmMinusOne; i++) {
       delete xSections[i];
     }
@@ -101 +103 @@
 
   G4int i;
-  G4int n = nElmMinusOne + 1;
-  G4double* xsec = new G4double[n];  
 
   // loop over bins
@@ -110 +110 @@
     cross = 0.0;
     //G4cout << "j= " << j << " e(MeV)= " << e/MeV << G4endl;
-    for (i=0; i<n; i++) {
+    for (i=0; i<=nElmMinusOne; i++) {
       cross += theAtomNumDensityVector[i]*      
         model->ComputeCrossSectionPerAtom(part, (*theElementVector)[i], e, 
                                           cutEnergy, e);
-      xsec[i] = cross;
-    }
-    if(DBL_MIN >= cross) cross = 1.0;
-    // normalise cross section sum 
-    for (i=0; i<nElmMinusOne; i++) {
-      xSections[i]->PutValue(j, xsec[i]/cross);
-      //G4cout << "i= " << i << " xs= " << xsec[i]/cross << G4endl;
+      xSections[i]->PutValue(j, cross);
     }
   }
-  delete [] xsec;
+
+  // xSections start from null, so use probabilities from the next bin
+  if(DBL_MIN >= (*xSections[nElmMinusOne])[0]) {
+    for (i=0; i<=nElmMinusOne; i++) {
+      xSections[i]->PutValue(0, (*xSections[i])[1]);
+    }
+  }
+  // xSections ends with null, so use probabilities from the previous bin
+  if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins-1]) {
+    for (i=0; i<=nElmMinusOne; i++) {
+      xSections[i]->PutValue(nbins-1, (*xSections[i])[nbins-2]);
+    }
+  }
+  // perform normalization
+  for(G4int j=0; j<nbins; j++) {
+    cross = (*xSections[nElmMinusOne])[j];
+    // only for positive X-section 
+    if(cross > DBL_MIN) {
+      for (i=0; i<nElmMinusOne; i++) {
+        xSections[i]->PutValue(j, (*xSections[i])[j]/cross);
+      }
+    }
+  }
 }
 
@@ -139 +155 @@
     }
   }  
-  G4cout << "Last Element in element vector" 
+  G4cout << "Last Element in element vector " 
          << (*theElementVector)[nElmMinusOne]->GetName() 
          << G4endl;