Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.cc,v 1.40 2007/11/09 11:35:54 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4EmModelManager.cc,v 1.46 2008/10/13 14:56:56 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7878//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    7979
    80 G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& list, G4DataVector& lowE)
     80G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx,
     81                               G4DataVector& lowE, const G4Region* reg)
    8182{
    8283  nModelsForRegion      = nMod;
    8384  theListOfModelIndexes = new G4int [nModelsForRegion];
    84   lowKineticEnergy      = new G4double [nModelsForRegion];
     85  lowKineticEnergy      = new G4double [nModelsForRegion+1];
    8586  for (G4int i=0; i<nModelsForRegion; i++) {
    86     theListOfModelIndexes[i] = list[i];
     87    theListOfModelIndexes[i] = indx[i];
    8788    lowKineticEnergy[i] = lowE[i];
    8889  }
     90  lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion];
     91  theRegion = reg;
    8992}
    9093
    … …  
    99102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    100103
    101 #include "G4LossTableManager.hh"
    102104#include "G4Step.hh"
    103105#include "G4ParticleDefinition.hh"
    … …  
    107109#include "G4MaterialCutsCouple.hh"
    108110#include "G4ProductionCutsTable.hh"
    109 #include "G4Region.hh"
    110111#include "G4RegionStore.hh"
    111112#include "G4Gamma.hh"
    112113#include "G4Positron.hh"
     114#include "G4UnitsTable.hh"
    113115
    114116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    128130  regions.clear();
    129131  orderOfModels.clear();
    130   upperEkin.clear();
    131132  maxCutInRange    = 12.*cm;
    132133  maxSubCutInRange = 0.7*mm;
    … …  
    139140G4EmModelManager::~G4EmModelManager()
    140141{
    141   G4int i,j;
    142142  Clear();
     143}
     144
     145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     146
     147void G4EmModelManager::Clear()
     148{
    143149  if(1 < verboseLevel) {
    144     G4cout << "G4EmModelManager:: delete models ";
    145     if(particle) G4cout << " for " << particle->GetParticleName();
    146     G4cout << " nModels=" << nEmModels <<G4endl;
    147   }
    148 
    149   for(i = 0; i<nEmModels; i++) {
    150     orderOfModels[i] = 1;
    151   }
    152   for(i = 0; i<nEmModels; i++) {
    153     if (orderOfModels[i]) {
    154       orderOfModels[i] = 0;
    155       for(j = i+1; j<nEmModels; j++) {
    156         if(models[i] == models[j]) orderOfModels[j] = 0;
    157       }
    158       G4String nam = models[i]->GetName();
    159       if(nam != "PAI" && nam != "PAIModel" ) delete models[i];
    160     }
    161   }
    162   for(i = 0; i<nEmModels; i++) {
    163     orderOfModels[i] = 1;
    164   }
    165   for(i = 0; i<nEmModels; i++) {
    166     if (orderOfModels[i]) {
    167       orderOfModels[i] = 0;
    168       for(j = i+1; j<nEmModels; j++) {
    169         if(flucModels[i] == flucModels[j]) orderOfModels[j] = 0;
    170       }
    171       delete flucModels[i];
    172     }
    173   }
    174   if(1 < verboseLevel)
    175     G4cout << "G4EmModelManager:: models are deleted!" << G4endl;
    176 }
    177 
    178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    179 
    180 void G4EmModelManager::Clear()
    181 {
    182   if(1 < verboseLevel) {
    183     G4cout << "G4EmModelManager::Clear()";
    184     if(particle) G4cout << " for " << particle->GetParticleName();
    185     G4cout << G4endl;
     150    G4cout << "G4EmModelManager::Clear()" << G4endl;
    186151  }
    187152
    188153  theCuts.clear();
    189154  theSubCuts.clear();
    190   upperEkin.clear();
    191155  if(idxOfRegionModels) delete [] idxOfRegionModels;
    192156  if(setOfRegionModels && nRegions) {
    … …  
    215179  orderOfModels.push_back(num);
    216180  p->DefineForRegion(r);
    217   if (nEmModels>0) {
    218     G4int idx = nEmModels;
    219     do {idx--;} while (idx && num < orderOfModels[idx]);
    220     if (num >= orderOfModels[idx] && num <= orderOfModels[idx+1]) idx++;
    221     if (idx < nEmModels) {
    222       models[nEmModels] = models[idx];
    223       flucModels[nEmModels] = flucModels[idx];
    224       regions[nEmModels] = regions[idx];
    225       orderOfModels[nEmModels] = orderOfModels[idx];
    226       models[idx] = p;
    227       flucModels[idx] = fm;
    228       regions[idx] = r;
    229       orderOfModels[idx] = num;
    230     }
    231   }
    232181  nEmModels++;
    233182}
    … …  
    254203//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    255204
    256 G4VEmModel* G4EmModelManager::GetModel(G4int i)
     205G4VEmModel* G4EmModelManager::GetModel(G4int i, G4bool ver)
    257206{
    258207  G4VEmModel* m = 0;
    259   if(i >= 0 && i < nEmModels) m = models[i];
    260   else if(verboseLevel > 0)
     208  if(i >= 0 && i < nEmModels) {m = models[i];}
     209  else if(verboseLevel > 0 && ver) {
    261210    G4cout << "G4EmModelManager::GetModel WARNING: "
    262211           << "index " << i << " is wrong Nmodels= "
    263            << nEmModels << G4endl;
     212           << nEmModels;
     213    if(particle) G4cout << " for " << particle->GetParticleName();
     214    G4cout<< G4endl;
     215  }
    264216  return m;
    265217}
    … …  
    269221const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
    270222                                                 const G4ParticleDefinition* sp,
    271                                                        G4double theMinSubRange,
    272                                                       G4int val)
     223                                                 G4double theMinSubRange,
     224                                                G4int val)
    273225{
    274226  verboseLevel = val;
    … …  
    293245  // Identify the list of regions with different set of models
    294246  nRegions = 1;
    295   std::vector<const G4Region*> set;
    296   set.push_back(world);
     247  std::vector<const G4Region*> setr;
     248  setr.push_back(world);
    297249  G4bool isWorld = false;
    298250
    … …  
    306258      if (nRegions>1) {
    307259        for (G4int j=1; j<nRegions; j++) {
    308           if ( r == set[j] ) newRegion = false;
     260          if ( r == setr[j] ) newRegion = false;
    309261        }
    310262      }
    311263      if (newRegion) {
    312         set.push_back(r);
     264        setr.push_back(r);
    313265        nRegions++;
    314266      }
    … …  
    322274  idxOfRegionModels[numOfCouples] = 0;
    323275  setOfRegionModels = new G4RegionModels*[nRegions];
    324   upperEkin.resize(nEmModels);
     276
     277  std::vector<G4int>    modelAtRegion(nEmModels);
     278  std::vector<G4int>    modelOrd(nEmModels);
     279  G4DataVector          eLow(nEmModels);
     280  G4DataVector          eHigh(nEmModels);
     281  G4int nmax = nEmModels;
    325282
    326283  // Order models for regions
    327284  for (G4int reg=0; reg<nRegions; reg++) {
    328     const G4Region* region = set[reg];
    329 
     285    const G4Region* region = setr[reg];
    330286    G4int n = 0;
    331 
    332     std::vector<G4int>    modelAtRegion;
    333     G4DataVector            eLow;
    334     G4DataVector            eHigh;
    335     modelAtRegion.clear();
    336     eLow.clear();
    337     eHigh.clear();
    338287
    339288    if(isWorld || 0 < reg) {
    … …  
    346295          G4double tmin = model->LowEnergyLimit();
    347296          G4double tmax = model->HighEnergyLimit();
    348           if (n>0) tmin = std::max(tmin, eHigh[n-1]);
     297          G4int  ord    = orderOfModels[ii];
     298          G4bool push   = true;
     299          G4bool insert = false;
     300          G4int  idx    = n;
    349301
    350302          if(1 < verboseLevel) {
    … …  
    355307                 << " tmin(MeV)= " << tmin/MeV
    356308                 << "; tmax(MeV)= " << tmax/MeV
    357                  << "; order= " << orderOfModels[ii]
     309                 << "; order= " << ord
    358310                 << G4endl;
    359311          }
    360 
    361           if (tmin < tmax) {
    362             modelAtRegion.push_back(ii);
    363             eLow.push_back(tmin);
    364             eHigh.push_back(tmax);
    365             upperEkin[ii] = tmax;
    366             n++;
     312       
     313          if (n == 0) n++;
     314          else {
     315            tmin = std::min(tmin, eHigh[n-1]);
     316            if(tmin >= tmax) push = false;
     317            else {
     318
     319              // high energy model
     320              if(tmin == eHigh[n-1] && tmax > eHigh[n-1]) n++;
     321              else if (tmax > eHigh[n-1]) {
     322
     323                // compare order of models
     324                for(G4int k = n-1; k>=0; k--) {
     325       
     326                  if (ord >= modelOrd[k]) {
     327                    tmin = std::max(tmin, eHigh[k]);           
     328                    if(k < n-1) n = k + 2;
     329                    break;
     330                  } else if (tmin > eLow[k]) {
     331                    eHigh[k] = tmin;
     332                    n = k + 2;
     333                    break;
     334                  } else if (tmin == eLow[k]) {
     335                    n = k + 1;
     336                    break;
     337                  }
     338                }
     339                if(tmin < eLow[0]) n = 1;
     340                idx = n - 1;
     341
     342                // low energy model
     343              } else {
     344
     345                tmax = std::max(tmax, eLow[0]);
     346                insert = true;
     347                push = false;
     348                idx = 0;
     349                if(tmax <= eLow[0]) tmax = eLow[0];
     350                else {
     351
     352                  for(G4int k=0; k<n; k++) {
     353       
     354                    if (ord >= modelOrd[k]) {
     355                      if(k == 0) {
     356                        if(tmin < eLow[0]) tmax = eLow[0];
     357                        else insert = false;
     358                        break;
     359                      } else {
     360                        insert = false;
     361                        break;
     362                      }
     363                    } else if(tmax < eHigh[k]) {
     364                      idx = k;
     365                      if(k > 0) tmin = eLow[k];
     366                      eLow[k] = tmax;               
     367                      break;
     368                    } else if(tmax == eHigh[k]) {           
     369                      insert = false;
     370                      push = true;
     371                      idx = k;
     372                      if(k > 0) tmin = eLow[k];
     373                      else tmin = std::min(tmin,eLow[0]);
     374                      break;
     375                    } else {
     376                      modelAtRegion[k] = ii;
     377                      modelOrd[k] = ord;
     378                      if(k == 0) eLow[idx] = std::min(tmin,eLow[0]);
     379                    }
     380                  }
     381                }
     382                if(insert && idx < n) n++;
     383                else insert = false;
     384              }
     385            }
     386          }
     387          if(n > nmax) {
     388            nmax = n;
     389            modelAtRegion.resize(nmax);
     390            modelOrd.resize(nmax);
     391            eLow.resize(nmax);
     392            eHigh.resize(nmax);
     393          }
     394          if(insert) {
     395            for(G4int k=n-2; k>=idx; k--) {           
     396              modelAtRegion[k+1] = modelAtRegion[k];
     397              modelOrd[k+1] = modelOrd[k];
     398              eLow[k+1]  = eLow[k];
     399              eHigh[k+1] = eHigh[k];
     400            }
     401          }
     402          if (push || insert) {
     403            modelAtRegion[idx] = ii;
     404            modelOrd[idx] = ord;
     405            eLow[idx]  = tmin;
     406            eHigh[idx] = tmax;
    367407          }
    368408        }
    … …  
    374414      eLow.push_back(0.0);
    375415      eHigh.push_back(DBL_MAX);
    376       upperEkin.push_back(DBL_MAX);
    377416    }
    378417    eLow[0] = 0.0;
     418    if(n >= nmax) eLow.resize(nmax+1);
     419    eLow[n] = eHigh[n-1];
    379420
    380421    if(1 < verboseLevel) {
    … …  
    382423      if (region) G4cout << region->GetName();
    383424      G4cout << ">  Elow(MeV)= ";
    384       for(G4int ii=0; ii<n; ii++) {G4cout << eLow[ii]/MeV << " ";}
     425      for(G4int ii=0; ii<=n; ii++) {G4cout << eLow[ii]/MeV << " ";}
    385426      G4cout << G4endl;
    386427    }
    387     G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow);
     428    G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
    388429    setOfRegionModels[reg] = rm;
    389430  }
    … …  
    399440 
    400441    G4int reg = nRegions;
    401     do {reg--;} while (reg>0 && pcuts != (set[reg]->GetProductionCuts()));
     442    do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
    402443    idxOfRegionModels[i] = reg;
    403444
    … …  
    474515{
    475516
    476   // vectors to provide continues dE/dx
    477   G4DataVector factor;
    478   G4DataVector dedxLow;
    479   G4DataVector dedxHigh;
    480 
    481517  G4double e;
    482518
    … …  
    491527    G4cout << "G4EmModelManager::FillDEDXVector() for "
    492528           << couple->GetMaterial()->GetName()
    493            << "  Ecut(MeV)= " << cut
    494            << "  Esubcut(MeV)= " << subcut
     529           << "  cut(MeV)= " << cut
     530           << "  subcut(MeV)= " << subcut
    495531           << "  Type " << tType
    496532           << "  for " << particle->GetParticleName()
    … …  
    501537  const G4RegionModels* regModels = setOfRegionModels[reg];
    502538  G4int nmod = regModels->NumberOfModels();
    503   factor.resize(nmod);
    504   dedxLow.resize(nmod);
    505   dedxHigh.resize(nmod);
     539
     540  // vectors to provide continues dE/dx
     541  G4DataVector factor(nmod);
     542  G4DataVector eLow(nmod+1);
     543  G4DataVector dedxLow(nmod);
     544  G4DataVector dedxHigh(nmod);
    506545
    507546  if(1 < verboseLevel) {
    … …  
    512551  }
    513552
    514 
    515553  // calculate factors to provide continuity of energy loss
     554
     555
    516556  factor[0] = 1.0;
    517557  G4int j;
    … …  
    519559  G4int totBinsLoss = aVector->GetVectorLength();
    520560
    521   dedxLow[0]  = 0.0;
    522 
    523   e = upperEkin[regModels->ModelIndex(0)];
     561  dedxLow[0] = 0.0;
     562  eLow[0]    = 0.0;
     563
     564  e = regModels->LowEdgeEnergy(1);
     565  eLow[1]    = e;
    524566  G4VEmModel* model = models[regModels->ModelIndex(0)];
    525567  dedxHigh[0] = 0.0;
    526568  if(model && cut > subcut) {
    527569    dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
    528     if(subcut > 0.0)
     570    if(subcut > 0.0) {
    529571      dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
     572    }
    530573  }
    531574  if(nmod > 1) {
    532575    for(j=1; j<nmod; j++) {
    533576
    534       e = upperEkin[regModels->ModelIndex(j-1)];
     577      e = regModels->LowEdgeEnergy(j);
     578      eLow[j]   = e;
    535579      G4int idx = regModels->ModelIndex(j);
    536580
    537581      dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
    538       if(subcut > 0.0)
     582      if(subcut > 0.0) {
    539583        dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
     584      }
    540585      if(subcut == cut) dedxLow[j] = 0.0;
    541586
    542       e = upperEkin[idx];
     587      e = regModels->LowEdgeEnergy(j+1);
     588      eLow[j+1] = e;
    543589      dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
    544       if(subcut > 0.0)
     590      if(subcut > 0.0) {
    545591        dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
     592      }
    546593      if(subcut == cut) dedxHigh[j] = 0.0;
    547594    }
     595    if(1 < verboseLevel) {
     596      G4cout << " model #0" 
     597             << "  dedx(" << eLow[0] << ")=  " << dedxLow[0]
     598             << "  dedx(" << eLow[1] << ")=  " << dedxHigh[0]
     599             << G4endl;
     600    }
    548601
    549602    for(j=1; j<nmod; j++) {
    550       if(dedxLow[j] > 0.0) factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0);
    551       else                 factor[j] = 0.0;
     603      if(dedxLow[j] > 0.0) {
     604        factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0)*eLow[j];
     605      } else  factor[j] = 0.0;
     606      if(1 < verboseLevel) {
     607        G4cout << " model #" << j
     608               << "  dedx(" << eLow[j] << ")=  " << dedxLow[j]
     609               << "  dedx(" << eLow[j+1] << ")=  " << dedxHigh[j]
     610               << "  factor= " << factor[j]/eLow[j]
     611               << G4endl;
     612      }
    552613    }
    553614
    … …  
    565626      // Choose a model of energy losses
    566627    G4int k = 0;
    567     if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
     628    if (nmod > 1 && e > eLow[1]) {
    568629      do {
    569630        k++;
    570         fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
    571       } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
     631        fac *= (1.0 + factor[k]/e);
     632      } while (k+1 < nmod && e > eLow[k+1]);
    572633    }
    573634
    … …  
    602663                                        G4EmTableType tType)
    603664{
    604   // vectors to provide continues cross section
    605   G4DataVector factor;
    606   G4DataVector sigmaLow;
    607   G4DataVector sigmaHigh;
    608 
    609665  G4double e;
    610666
    … …  
    630686  const G4RegionModels* regModels = setOfRegionModels[reg];
    631687  G4int nmod = regModels->NumberOfModels();
    632   factor.resize(nmod);
    633   sigmaLow.resize(nmod);
    634   sigmaHigh.resize(nmod);
     688
     689  // vectors to provide continues dE/dx
     690  G4DataVector factor(nmod);
     691  G4DataVector eLow(nmod+1);
     692  G4DataVector sigmaLow(nmod);
     693  G4DataVector sigmaHigh(nmod);
    635694
    636695  if(2 < verboseLevel) {
    … …  
    644703  G4int totBinsLambda = aVector->GetVectorLength();
    645704
    646   sigmaLow[0]  = 0.0;
    647 
    648   e = upperEkin[regModels->ModelIndex(0)];
     705  sigmaLow[0] = 0.0;
     706  eLow[0]     = 0.0;
     707
     708  e = regModels->LowEdgeEnergy(1);
     709  eLow[1]     = e;
    649710  G4VEmModel* model = models[regModels->ModelIndex(0)];
    650711  sigmaHigh[0] = 0.0;
    … …  
    668729    for(j=1; j<nmod; j++) {
    669730
    670       e  = upperEkin[regModels->ModelIndex(j-1)];
    671       sigmaLow[j] =
    672         models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
    673       e  = upperEkin[regModels->ModelIndex(j)];
    674       sigmaHigh[j] =
    675         models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
     731      e  = regModels->LowEdgeEnergy(j);
     732      eLow[j]   = e;
     733      G4int idx = regModels->ModelIndex(j);
     734
     735      sigmaLow[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
     736      e  = regModels->LowEdgeEnergy(j+1);
     737      eLow[j+1] = e;
     738      sigmaHigh[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
     739    }
     740    if(1 < verboseLevel) {
     741      G4cout << " model #0" 
     742             << "  sigma(" << eLow[0] << ")=  " << sigmaLow[0]
     743             << "  sigma(" << eLow[1] << ")=  " << sigmaHigh[0]
     744             << G4endl;
     745    }
     746    for(j=1; j<nmod; j++) {
     747      if(sigmaLow[j] > 0.0) {
     748        factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0)*eLow[j];
     749      } else  factor[j] = 0.0;
    676750      if(1 < verboseLevel) {
    677         G4cout << " model #" << j << "   eUp= " << e
    678                << "  sigmaUp= " << sigmaHigh[j]
    679                << "  sigmaDown= " << sigmaLow[j]
     751        G4cout << " model #" << j
     752               << "  sigma(" << eLow[j] << ")=  " << sigmaLow[j]
     753               << "  sigma(" << eLow[j+1] << ")=  " << sigmaHigh[j]
     754               << "  factor= " << factor[j]/eLow[j]
    680755               << G4endl;
    681756      }
    682     }
    683     for(j=1; j<nmod; j++) {
    684       if(sigmaLow[j] > 0.0) factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0);
    685       else                  factor[j] = 0.0;
    686757    }
    687758  }
    … …  
    695766    G4int k = 0;
    696767    G4double fac = 1.0;
    697     if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
     768    if (nmod > 1 && e > eLow[1]) {
    698769      do {
    699770        k++;
    700         fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
    701       } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
     771        fac *= (1.0 + factor[k]/e);
     772      } while ( k+1 < nmod && e > eLow[k+1] );
    702773    }
    703774
    … …  
    721792
    722793//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     794
     795void G4EmModelManager::DumpModelList(G4int verb)
     796{
     797  if(verb == 0) return;
     798  for(G4int i=0; i<nRegions; i++) {
     799    G4RegionModels* r = setOfRegionModels[i];
     800    const G4Region* reg = r->Region();
     801    if(verb > 1 || nRegions > 1) {
     802    }
     803    G4int n = r->NumberOfModels(); 
     804    if(verb > 1 || n > 0) {
     805      G4cout << "      ===== EM models for the G4Region  " << reg->GetName()
     806             << " ======" << G4endl;;
     807      for(G4int j=0; j<n; j++) {
     808        const G4VEmModel* m = models[r->ModelIndex(j)];
     809        G4cout << std::setw(20);
     810        G4cout << m->GetName() << " :     Emin= "
     811               << std::setw(10) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
     812               << "        Emax=   "
     813               << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy")
     814               << G4endl;
     815      } 
     816    }
     817  }
     818}
     819
     820//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracChangeset for help on using the changeset viewer.