Ignore:
Timestamp:
May 28, 2009, 4:26:57 PM (15 years ago)
Author:
garnier
Message:

maj sur la beta de geant 4.9.3

Location:
trunk/source/processes/electromagnetic/lowenergy/src
Files:
60 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/lowenergy/src/G4AtomicDeexcitation.cc

    r1007 r1055  
    2626//
    2727// $Id: G4AtomicDeexcitation.cc,v 1.11
    28 // GEANT4 tag $Name: geant4-09-02 $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Authors: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
    … …  
    111111    }
    112112 
    113   // Look this in a particular way: only one auger emitted! //
     113  // Look this in a particular way: only one auger emitted! // ????
    114114  while (provShellId > -2);
    115115 
    … …  
    297297//            G4cout << "G4AtomicDeexcitation warning: No Auger transition found" <<  G4endl;
    298298//            G4cout << "Absorbed enrgy deposited locally" << G4endl;
    299               return 0;
    300 //            //  G4Exception("G4AtomicDeexcitation: No Auger transition found");
     299              G4Exception("G4AtomicDeexcitation: No Auger transition found");
     300              return 0;
    301301            }
    302302        }
    … …  
    394394      // G4int augerOriginatingShellId = 0;
    395395     
    396       G4int numberOfPossibleAuger =
    397           (anAugerTransition->AugerTransitionProbabilities(transitionRandomShellId))->size();
     396      G4int numberOfPossibleAuger = 0;
     397      numberOfPossibleAuger = anAugerTransition->AugerTransitionProbabilities(transitionRandomShellId)->size();
     398
     399
    398400      G4bool foundFlag = false;
    399401

  • trunk/source/processes/electromagnetic/lowenergy/src/G4AugerData.cc

    r961 r1055  
    486486      //        G4cout << "G4AugerData for Element no. " << element << " are loaded" << G4endl;
    487487      // G4cout << "G4AugerData for Element no. " << element << " are loaded" << G4endl;
    488       G4cout << "AugerTransitionTable complete"<< G4endl;
     488      //G4cout << "AugerTransitionTable complete"<< G4endl;
    489489    }
    490490}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4BremsstrahlungCrossSectionHandler.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BremsstrahlungCrossSectionHandler.cc,v 1.9 2006/06/29 19:38:42 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4BremsstrahlungCrossSectionHandler.cc,v 1.10 2009/03/03 11:19:18 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4141// 10.10.2001 MGP Revision to improve code quality and consistency with design
    4242// 21.01.2003 VI  cut per region
     43// 03.03.2009 LP  Added public method to make a easier migration of
     44//                G4LowEnergyBremsstrahlung to G4LivermoreBremsstrahlungModel
    4345//
    4446// -------------------------------------------------------------------
    … …  
    5557#include "G4ProductionCutsTable.hh"
    5658
     59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     60
    5761G4BremsstrahlungCrossSectionHandler::G4BremsstrahlungCrossSectionHandler(const G4VEnergySpectrum* spec,
    5862                                                                         G4VDataSetAlgorithm* )
    … …  
    6266}
    6367
     68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    6469
    6570G4BremsstrahlungCrossSectionHandler::~G4BremsstrahlungCrossSectionHandler()
    … …  
    6873}
    6974
     75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    7076
    7177std::vector<G4VEMDataSet*>*
    … …  
    127133  return set;
    128134}
     135
     136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     137
     138G4double G4BremsstrahlungCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
     139                                                                                      G4double cutEnergy,
     140                                                                                      G4int Z)
     141{
     142  G4double value = 0.;
     143  if(energy > cutEnergy)
     144    {
     145      G4double elemCs = FindValue(Z, energy);
     146      value  = theBR->Probability(Z,cutEnergy, energy, energy);   
     147      value *= elemCs;
     148    }
     149  return value;
     150}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionChargeDecrease.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionChargeDecrease.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionChargeDecrease.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionChargeDecrease.hh"
    … …  
    7979  }
    8080
     81   G4cout << G4endl;
     82   G4cout << "*******************************************************************************" << G4endl;
     83   G4cout << "*******************************************************************************" << G4endl;
     84   G4cout << "   The class G4CrossSectionChargeDecrease is NOT SUPPORTED ANYMORE. " << G4endl;
     85   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     86   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     87   G4cout << "*******************************************************************************" << G4endl;
     88   G4cout << "*******************************************************************************" << G4endl;
     89   G4cout << G4endl;
    8190}
    8291

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionChargeIncrease.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionChargeIncrease.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionChargeIncrease.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionChargeIncrease.hh"
    … …  
    7979  }
    8080
     81   G4cout << G4endl;
     82   G4cout << "*******************************************************************************" << G4endl;
     83   G4cout << "*******************************************************************************" << G4endl;
     84   G4cout << "   The class G4CrossSectionChargeIncrease is NOT SUPPORTED ANYMORE. " << G4endl;
     85   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     86   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     87   G4cout << "*******************************************************************************" << G4endl;
     88   G4cout << "*******************************************************************************" << G4endl;
     89   G4cout << G4endl;
    8190}
    8291

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionElasticChampion.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionElasticChampion.cc,v 1.4 2008/12/05 11:58:16 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionElasticChampion.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828// -------------------------------------------------------------------
    2929
    … …  
    5757    G4Exception("G4CrossSectionElasticChampion constructor: electron is not defined");
    5858  }
     59
     60   G4cout << G4endl;
     61   G4cout << "*******************************************************************************" << G4endl;
     62   G4cout << "*******************************************************************************" << G4endl;
     63   G4cout << "   The class G4CrossSectionElasticChampion is NOT SUPPORTED ANYMORE. " << G4endl;
     64   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     65   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     66   G4cout << "*******************************************************************************" << G4endl;
     67   G4cout << "*******************************************************************************" << G4endl;
     68   G4cout << G4endl;
    5969}
    6070

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionElasticScreenedRutherfordHE.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionElasticScreenedRutherfordHE.cc,v 1.2 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionElasticScreenedRutherfordHE.cc,v 1.3 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionElasticScreenedRutherfordHE.hh"
    … …  
    3535  lowEnergyLimit = 200. * eV;
    3636  highEnergyLimit = 10. * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4CrossSectionElasticScreenedRutherfordHE is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionElasticScreenedRutherfordLE.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionElasticScreenedRutherfordLE.cc,v 1.2 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionElasticScreenedRutherfordLE.cc,v 1.3 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionElasticScreenedRutherfordLE.hh"
    … …  
    3535  lowEnergyLimit = 0 * eV;
    3636  highEnergyLimit = 200 * eV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4CrossSectionElasticScreenedRutherfordLE is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionExcitationBorn.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionExcitationBorn.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionExcitationBorn.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionExcitationBorn.hh"
    … …  
    3636  highEnergyLimit = 10 * MeV;
    3737  table = 0;
     38
     39   G4cout << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "*******************************************************************************" << G4endl;
     42   G4cout << "   The class G4CrossSectionExcitationBorn is NOT SUPPORTED ANYMORE. " << G4endl;
     43   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     44   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << "*******************************************************************************" << G4endl;
     47   G4cout << G4endl;
    3848}
    3949

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionExcitationEmfietzoglou.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionExcitationEmfietzoglou.cc,v 1.5 2008/12/05 11:58:16 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionExcitationEmfietzoglou.cc,v 1.6 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionExcitationEmfietzoglou.hh"
    … …  
    3535  lowEnergyLimit = 8.23 * eV;
    3636  highEnergyLimit = 10. * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4CrossSectionExcitationEmfietzoglou is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionExcitationMillerGreen.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionExcitationMillerGreen.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionExcitationMillerGreen.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionExcitationMillerGreen.hh"
    … …  
    9292  }
    9393   
     94
     95   G4cout << G4endl;
     96   G4cout << "*******************************************************************************" << G4endl;
     97   G4cout << "*******************************************************************************" << G4endl;
     98   G4cout << "   The class G4CrossSectionExcitationMillerGreen is NOT SUPPORTED ANYMORE. " << G4endl;
     99   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     100   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     101   G4cout << "*******************************************************************************" << G4endl;
     102   G4cout << "*******************************************************************************" << G4endl;
     103   G4cout << G4endl;
    94104}
    95105

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionExcitationMillerGreenPartial.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionExcitationMillerGreenPartial.cc,v 1.2 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionExcitationMillerGreenPartial.cc,v 1.3 2009/01/20 07:40:53 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionExcitationMillerGreenPartial.hh"
    … …  
    115115  tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
    116116
     117  // SI - added protection
     118  if (tCorrected < waterExcitation.ExcitationEnergy(excitationLevel)) return 0;
     119  //
     120 
    117121  G4int z = 10;
    118122
    … …  
    158162  // ELECTRON CORRECTION
    159163 
    160   if ( particle == instance->GetIon("alpha++"))
     164  if ( particle == instance->GetIon("alpha++")||
     165       particle == G4Proton::ProtonDefinition())
     166       
    161167  {  while (i > 0)
    162168     {

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionIonisationBorn.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionIonisationBorn.cc,v 1.4 2008/08/20 14:51:48 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionIonisationBorn.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionIonisationBorn.hh"
    … …  
    8282    G4Exception("G4CrossSectionIonisationBorn Constructor: proton is not defined");
    8383  }
     84
     85   G4cout << G4endl;
     86   G4cout << "*******************************************************************************" << G4endl;
     87   G4cout << "*******************************************************************************" << G4endl;
     88   G4cout << "   The class G4CrossSectionIonisationBorn is NOT SUPPORTED ANYMORE. " << G4endl;
     89   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     90   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     91   G4cout << "*******************************************************************************" << G4endl;
     92   G4cout << "*******************************************************************************" << G4endl;
     93   G4cout << G4endl;
    8494}
    8595

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionIonisationBornPartial.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionIonisationBornPartial.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionIonisationBornPartial.cc,v 1.5 2009/01/20 07:40:53 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionIonisationBornPartial.hh"
    … …  
    3333G4CrossSectionIonisationBornPartial::G4CrossSectionIonisationBornPartial()
    3434{
    35   lowEnergyLimitDefault = 25 * eV;
     35  lowEnergyLimitDefault = 12.61 * eV;
    3636  highEnergyLimitDefault = 30 * keV;
    3737
    … …  
    5252    tableFile[electron] = fileElectron;
    5353
    54     lowEnergyLimit[electron] = 25. * eV;
     54    lowEnergyLimit[electron] = 12.61 * eV;
    5555    highEnergyLimit[electron] = 30. * keV;
    5656

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionIonisationRudd.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CrossSectionIonisationRudd.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4CrossSectionIonisationRudd.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4CrossSectionIonisationRudd.hh"
    … …  
    154154    G4Exception("G4CrossSectionIonisationRudd Constructor: helium is not defined");
    155155  }
     156
     157   G4cout << G4endl;
     158   G4cout << "*******************************************************************************" << G4endl;
     159   G4cout << "*******************************************************************************" << G4endl;
     160   G4cout << "   The class G4CrossSectionIonisationRudd is NOT SUPPORTED ANYMORE. " << G4endl;
     161   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     162   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     163   G4cout << "*******************************************************************************" << G4endl;
     164   G4cout << "*******************************************************************************" << G4endl;
     165   G4cout << G4endl;
     166
    156167}
    157168

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateChargeDecrease.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateChargeDecrease.cc,v 1.3 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateChargeDecrease.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateChargeDecrease.hh"
    … …  
    3535  lowEnergyLimit = 1 * keV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4FinalStateChargeDecrease is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848
    … …  
    7080  }
    7181 
     82  //SI - Added protection against total energy deposit
     83  product.DoNotDepositEnergy();
     84  //
    7285  product.KillPrimaryParticle();
     86
    7387  product.AddEnergyDeposit(waterBindingEnergy);
    74  
     88
    7589  G4DynamicParticle* aSecondary = new G4DynamicParticle(OutgoingParticleDefinition(definition, finalStateIndex),
    7690                                                        track.GetDynamicParticle()->GetMomentumDirection(),

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateChargeIncrease.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateChargeIncrease.cc,v 1.3 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateChargeIncrease.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateChargeIncrease.hh"
    … …  
    3535  lowEnergyLimit = 1 * keV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class 4FinalStateChargeIncrease is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848
    … …  
    4858  product.Clear();
    4959
     60  //SI - Added protection against total energy deposit
     61  product.DoNotDepositEnergy();
     62  //
    5063  product.KillPrimaryParticle();
    5164  product.AddEnergyDeposit(0.);
     65
    5266  G4ParticleDefinition* definition = track.GetDefinition();
    5367 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateElasticBrennerZaider.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateElasticBrennerZaider.cc,v 1.8 2008/12/05 11:58:16 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateElasticBrennerZaider.cc,v 1.9 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateElasticBrennerZaider.hh"
    … …  
    6464  gamma100_200Coeff.push_back(-2.96264E-5);
    6565  gamma100_200Coeff.push_back(-1.20655E-7);
     66
     67   G4cout << G4endl;
     68   G4cout << "*******************************************************************************" << G4endl;
     69   G4cout << "*******************************************************************************" << G4endl;
     70   G4cout << "   The class G4FinalStateElasticBrennerZaider is NOT SUPPORTED ANYMORE. " << G4endl;
     71   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     72   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     73   G4cout << "*******************************************************************************" << G4endl;
     74   G4cout << "*******************************************************************************" << G4endl;
     75   G4cout << G4endl;
    6676}
    6777

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateElasticChampion.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateElasticChampion.cc,v 1.7 2008/12/10 18:25:28 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateElasticChampion.cc,v 1.8 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828// -------------------------------------------------------------------
    2929
    … …  
    8484      G4Exception("G4FinalStateElastiChampion : constructor: electron is not defined");
    8585  }
     86
     87   G4cout << G4endl;
     88   G4cout << "*******************************************************************************" << G4endl;
     89   G4cout << "*******************************************************************************" << G4endl;
     90   G4cout << "   The class G4FinalStateElastiChampion is NOT SUPPORTED ANYMORE. " << G4endl;
     91   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     92   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     93   G4cout << "*******************************************************************************" << G4endl;
     94   G4cout << "*******************************************************************************" << G4endl;
     95   G4cout << G4endl;
    8696}
    8797

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateElasticScreenedRutherford.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateElasticScreenedRutherford.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateElasticScreenedRutherford.cc,v 1.5 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateElasticScreenedRutherford.hh"
    … …  
    3535  lowEnergyLimit = 200 * eV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4FinalStateElasticScreenedRutherford is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateExcitationBorn.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateExcitationBorn.cc,v 1.3 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateExcitationBorn.cc,v 1.4 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateExcitationBorn.hh"
    … …  
    3535  lowEnergyLimit = 500 * keV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4FinalStateExcitationBorn is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateExcitationEmfietzoglou.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateExcitationEmfietzoglou.cc,v 1.5 2008/12/05 11:58:16 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateExcitationEmfietzoglou.cc,v 1.6 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateExcitationEmfietzoglou.hh"
    … …  
    3535  lowEnergyLimit = 8.23 * eV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4FinalStateExcitationEmfietzoglou is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateExcitationMillerGreen.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateExcitationMillerGreen.cc,v 1.3 2008/07/14 20:47:34 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateExcitationMillerGreen.cc,v 1.4 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateExcitationMillerGreen.hh"
    … …  
    3535  lowEnergyLimit = 10 * eV;
    3636  highEnergyLimit = 10 * MeV;
     37
     38   G4cout << G4endl;
     39   G4cout << "*******************************************************************************" << G4endl;
     40   G4cout << "*******************************************************************************" << G4endl;
     41   G4cout << "   The class G4FinalStateExcitationMillerGreen is NOT SUPPORTED ANYMORE. " << G4endl;
     42   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     43   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     44   G4cout << "*******************************************************************************" << G4endl;
     45   G4cout << "*******************************************************************************" << G4endl;
     46   G4cout << G4endl;
    3747}
    3848

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationBorn.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateIonisationBorn.cc,v 1.16 2008/12/06 13:47:12 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateIonisationBorn.cc,v 1.17 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateIonisationBorn.hh"
    … …  
    124124    G4Exception("G4FinalStateIonisationBorn Constructor: proton is not defined");
    125125  }
     126
     127   G4cout << G4endl;
     128   G4cout << "*******************************************************************************" << G4endl;
     129   G4cout << "*******************************************************************************" << G4endl;
     130   G4cout << "   The class G4FinalStateIonisationBorn is NOT SUPPORTED ANYMORE. " << G4endl;
     131   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     132   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     133   G4cout << "*******************************************************************************" << G4endl;
     134   G4cout << "*******************************************************************************" << G4endl;
     135   G4cout << G4endl;
    126136}
    127137

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4FinalStateIonisationRudd.cc,v 1.8 2008/08/20 14:51:48 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FinalStateIonisationRudd.cc,v 1.9 2009/05/02 15:07:47 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828
    2929#include "G4FinalStateIonisationRudd.hh"
    … …  
    7070  lowEnergyLimit[helium] = 1. * keV;
    7171  highEnergyLimit[helium] = 10. * MeV;
     72
     73   G4cout << G4endl;
     74   G4cout << "*******************************************************************************" << G4endl;
     75   G4cout << "*******************************************************************************" << G4endl;
     76   G4cout << "   The class G4FinalStateIonisationRudd is NOT SUPPORTED ANYMORE. " << G4endl;
     77   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     78   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     79   G4cout << "*******************************************************************************" << G4endl;
     80   G4cout << "*******************************************************************************" << G4endl;
     81   G4cout << G4endl;
    7282}
    7383

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateProduct.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4FinalStateProduct.cc,v 1.5 2007/11/09 20:11:04 pia Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4FinalStateProduct.cc,v 1.6 2009/01/20 07:50:28 sincerti Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
    … …  
    5353#include "G4ThreeVector.hh"
    5454
    55 G4FinalStateProduct::G4FinalStateProduct() : killStatus(false), isModified(false), localEnergyDeposit(0.), modifiedEnergy(0)
     55G4FinalStateProduct::G4FinalStateProduct() : killStatus(false), doNotDepositStatus(false), isModified(false), localEnergyDeposit(0.), modifiedEnergy(0)
    5656{
    5757  // empty
    … …  
    6868  // Reset object status
    6969  killStatus = false;
     70  doNotDepositStatus = false;
    7071  isModified = false;
    7172  localEnergyDeposit = 0.;
    … …  
    9899}
    99100 
     101void G4FinalStateProduct::DoNotDepositEnergy()
     102{
     103  doNotDepositStatus = true;
     104}
     105
    100106void G4FinalStateProduct::KillPrimaryParticle()
    101107{
     108 
    102109  // ---- MGP ---- To be added: Handle local energy deposit here
    103110  killStatus = true;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4IonParametrisedLossModel.cc

    r1007 r1055  
    3737// First implementation: 10. 11. 2008
    3838//
    39 // Modifications:
    40 //
     39// Modifications: 03. 02. 2009 - Bug fix iterators (AL)
     40//                11. 03. 2009 - Introduced new table handler (G4IonDEDXHandler)
     41//                               and modified method to add/remove tables
     42//                               (tables are now built in initialisation phase),
     43//                               Minor bug fix in ComputeDEDXPerVolume (AL)
     44//                11. 05. 2009 - Introduced scaling algorithm for heavier ions:
     45//                               G4IonDEDXScalingICRU73 (AL)
    4146//
    4247// Class description:
    … …  
    5560#include "G4MaterialStoppingICRU73.hh"
    5661#include "G4SimpleMaterialStoppingICRU73.hh"
     62#include "G4IronStoppingICRU73.hh"
     63#include "G4VIonDEDXTable.hh"
     64#include "G4VIonDEDXScalingAlgorithm.hh"
     65#include "G4IonDEDXScalingICRU73.hh"
    5766#include "G4BraggIonModel.hh"
    5867#include "G4BetheBlochModel.hh"
     68#include "G4ProductionCutsTable.hh"
    5969#include "G4ParticleChangeForLoss.hh"
    6070#include "G4LossTableManager.hh"
    … …  
    6373#include "Randomize.hh"
    6474
     75//#define PRINT_TABLE_BUILT
     76
     77
     78// #########################################################################
    6579
    6680G4IonParametrisedLossModel::G4IonParametrisedLossModel(
    … …  
    89103  betheBlochModel = new G4BetheBlochModel();
    90104
    91   // By default ICRU 73 stopping power tables are loaded
    92   AddDEDXTable<G4SimpleMaterialStoppingICRU73>();
    93   AddDEDXTable<G4MaterialStoppingICRU73>();
     105  // By default ICRU 73 stopping power tables are loaded:
     106
     107  // Ions with Z above between 19 and 21: Ar-40 data is used as basis
     108  // for stopping power scaling
     109  G4int ionZMin = 19;
     110  G4int ionZMax = 21;
     111  G4int refIonZ = 18;
     112  G4int refIonA = 40;
     113
     114  AddDEDXTable("ICRU73-elemmat",
     115           new G4SimpleMaterialStoppingICRU73,
     116           new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA));
     117
     118  // Ions with Z above 21: Fe-56 data is used as basis for stopping power
     119  // scaling
     120  ionZMin = 22;
     121  ionZMax = 102;
     122  refIonZ = 26;
     123  refIonA = 56;
     124
     125  AddDEDXTable("ICRU73-ironions",
     126               new G4IronStoppingICRU73,
     127               new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA));
     128
     129  // Compound materials: Ar-40 data is used as basis for stopping power
     130  // scaling (except for iron ions)
     131  ionZMin = 19;
     132  ionZMax = 102;
     133  refIonZ = 18;
     134  refIonA = 40;
     135
     136  G4IonDEDXScalingICRU73* scaling =
     137            new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA);
     138
     139  G4int ironIonAtomicNumber = 26;
     140  scaling -> AddException(ironIonAtomicNumber);
     141
     142  AddDEDXTable("ICRU73-compmat",
     143               new G4MaterialStoppingICRU73,
     144               scaling);
    94145
    95146  // The boundaries for the range tables are set
    … …  
    107158  dedxCacheMaterial = 0;
    108159  dedxCacheEnergyCut = 0;
    109   dedxCacheIter = lossTableList.begin();
     160  dedxCacheIter = lossTableList.end();
    110161  dedxCacheTransitionEnergy = 0.0; 
    111162  dedxCacheTransitionFactor = 0.0;
    … …  
    113164}
    114165
     166// #########################################################################
    115167
    116168G4IonParametrisedLossModel::~G4IonParametrisedLossModel() {
    … …  
    142194}
    143195
     196// #########################################################################
    144197
    145198G4double G4IonParametrisedLossModel::MinEnergyCut(
    … …  
    151204}
    152205
     206// #########################################################################
    153207
    154208void G4IonParametrisedLossModel::Initialise(
    … …  
    166220  dedxCacheMaterial = 0;
    167221  dedxCacheEnergyCut = 0;
    168   dedxCacheIter = lossTableList.begin();
     222  dedxCacheIter = lossTableList.end();
    169223  dedxCacheTransitionEnergy = 0.0; 
    170224  dedxCacheTransitionFactor = 0.0;
    … …  
    196250  cutEnergies.clear();
    197251  for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
     252
     253  // All dE/dx vectors are built
     254  const G4ProductionCutsTable* coupleTable=
     255                     G4ProductionCutsTable::GetProductionCutsTable();
     256  size_t nmbCouples = coupleTable -> GetTableSize();
     257
     258#ifdef PRINT_TABLE_BUILT
     259    G4cout << "G4IonParametrisedLossModel::Initialise():"
     260           << " Building dE/dx vectors:"
     261           << G4endl;         
     262#endif
     263
     264  for (size_t i = 0; i < nmbCouples; i++) {
     265
     266    const G4MaterialCutsCouple* couple =
     267                                     coupleTable -> GetMaterialCutsCouple(i);
     268
     269    const G4Material* material = couple -> GetMaterial();
     270    //    G4ProductionCuts* productionCuts = couple -> GetProductionCuts();
     271
     272    for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
     273
     274       LossTableList::iterator iter = lossTableList.begin();
     275       LossTableList::iterator iter_end = lossTableList.end();
     276
     277       for(;iter != iter_end; iter++) {
     278
     279          if(*iter == 0) {
     280              G4cout << "G4IonParametrisedLossModel::Initialise():"
     281                     << " Skipping illegal table."
     282                     << G4endl;         
     283          }
     284
     285          G4bool isApplicable =
     286                    (*iter) -> BuildDEDXTable(atomicNumberIon, material);
     287          if(isApplicable) {
     288
     289#ifdef PRINT_TABLE_BUILT
     290             G4cout << "  Atomic Number Ion = " << atomicNumberIon
     291                    << ", Material = " << material -> GetName()
     292                    << ", Table = " << (*iter) -> GetName()
     293                    << G4endl;     
     294#endif
     295             break;
     296          }
     297       }
     298    }
     299  }
    198300
    199301  // The particle change object is cast to G4ParticleChangeForLoss
    … …  
    221323}
    222324
     325// #########################################################################
    223326
    224327G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom(
    … …  
    285388}
    286389
     390// #########################################################################
    287391
    288392G4double G4IonParametrisedLossModel::CrossSectionPerVolume(
    … …  
    303407}
    304408
     409// #########################################################################
    305410
    306411G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume(
    … …  
    342447  LossTableList::iterator iter = dedxCacheIter;
    343448
    344   if(iter != lossTableList.begin()) {
     449  if(iter != lossTableList.end()) {
    345450
    346451     G4double transitionEnergy = dedxCacheTransitionEnergy;
    … …  
    369474        if(scaledTransitionEnergy >= lowEnergyLimit) {
    370475
    371            G4double factor = 1.0 + dedxCacheTransitionFactor /
    372                                                            kineticEnergy;
    373 
    374476           dEdx = betheBlochModel -> ComputeDEDXPerVolume(
    375477                                      material, genericIon,
    376478                                      scaledKineticEnergy, cutEnergy);
     479       
     480           dEdx *= chargeSquare;
     481
     482           dEdx += corrections -> ComputeIonCorrections(particle,
     483                                                 material, kineticEnergy);
     484
     485           G4double factor = 1.0 + dedxCacheTransitionFactor /
     486                                                           kineticEnergy;
     487
    377488           dEdx *= factor;
    378 
    379         }
    380         dEdx *= chargeSquare;
    381 
    382         dEdx += corrections -> ComputeIonCorrections(particle,
    383                                                  material, kineticEnergy);
     489        }
    384490     }
    385491  }
    … …  
    446552  if (dEdx < 0.0) dEdx = 0.0;
    447553
    448 #ifdef PRINT_DEBUG
    449 
    450   G4cout << "########################################################"
    451          << G4endl
    452          << "# G4IonParametrisedLossModel::ComputeDEDXPerVolume"
    453          << G4endl
    454          << "# Material =" << material -> GetName()
    455          << G4endl
    456          << "# Particle = " << particle -> GetParticleName()           
    457          << G4endl;
    458          << "# Cut energy (MeV) = " << cutEnergy/MeV 
    459          << G4endl;
    460  
    461   G4cout << "#"
    462          << std::setw(13) << std::right << "E(MeV)"
    463          << std::setw(14) << "dE/dx(keV/um)"
    464          << std::setw(14) << "d:dE/dx(keV/um)"
    465          << std::setw(14) << "(d:dE/dx)/dE/dx"
    466          << G4endl
    467          << "# ------------------------------------------------------"
    468          << G4endl;
    469 
    470   G4cout << std::setw(14) << std::right << kineticEnergy / MeV
    471          << std::setw(14) << (dEdx + dEdXDeltaRays) / keV * um
    472          << std::setw(14) << dEdXDeltaRays / keV * um
    473          << std::setw(14) << dEdXDeltaRays / (dEdx + dEdXDeltaRays) * 100.0
    474          << G4endl;
    475 #endif
    476 
    477554  return dEdx;
    478555}
    479556
     557// #########################################################################
    480558
    481559void G4IonParametrisedLossModel::PrintDEDXTable(
    … …  
    508586         << std::setw(13) << std::right << "(MeV)"
    509587         << std::setw(14) << "(MeV)"
    510          << std::setw(14) << "(MeV/mm)"
     588         << std::setw(14) << "(MeV/cm)"
    511589         << std::setw(14) << "(MeV*cm2/mg)"
    512590         << G4endl
    … …  
    535613      G4cout << std::setw(14) << std::right << energy / MeV
    536614             << std::setw(14) << energy / atomicMassNumber / MeV
    537              << std::setw(14) << dedx / MeV * mm
     615             << std::setw(14) << dedx / MeV * cm
    538616             << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
    539617             << G4endl;
    … …  
    541619}
    542620
     621// #########################################################################
     622
     623void G4IonParametrisedLossModel::PrintDEDXTableHandlers(
     624                   const G4ParticleDefinition* particle,  // Projectile (ion)
     625                   const G4Material* material,  // Absorber material
     626                   G4double lowerBoundary,      // Minimum energy per nucleon
     627                   G4double upperBoundary,      // Maximum energy per nucleon
     628                   G4int nmbBins,               // Number of bins
     629                   G4bool logScaleEnergy) {     // Logarithmic scaling of energy
     630
     631  LossTableList::iterator iter = lossTableList.begin();
     632  LossTableList::iterator iter_end = lossTableList.end();
     633
     634  for(;iter != iter_end; iter++) {
     635      G4bool isApplicable =  (*iter) -> IsApplicable(particle, material);
     636      if(isApplicable) { 
     637        (*iter) -> PrintDEDXTable(particle, material,
     638                                  lowerBoundary, upperBoundary,
     639                                  nmbBins,logScaleEnergy);
     640        break;
     641      }
     642  }
     643}
     644
     645// #########################################################################
    543646
    544647void G4IonParametrisedLossModel::SampleSecondaries(
    … …  
    641744}
    642745
     746// #########################################################################
    643747
    644748void G4IonParametrisedLossModel::UpdateDEDXCache(
    … …  
    672776
    673777     // If any table is applicable, the transition factor is computed:
    674      if(iter != lossTableList.begin()) {
     778     if(iter != lossTableList.end()) {
    675779
    676780        // Retrieving the transition energy from the parameterisation table
    … …  
    740844}
    741845
     846// #########################################################################
    742847
    743848void G4IonParametrisedLossModel::CorrectionsAlongStep(
    … …  
    767872  G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
    768873
     874  if(kineticEnergy == eloss) { return; }
     875
    769876  G4double cutEnergy = DBL_MAX;
    770877  size_t cutIndex = couple -> GetIndex();
    … …  
    777884  // If parameterization for ions is available the electronic energy loss
    778885  // is overwritten
    779   if(iter != lossTableList.begin()) {
     886  if(iter != lossTableList.end()) {
    780887
    781888     // The energy loss is calculated using the ComputeDEDXPerVolume function
    … …  
    832939
    833940     }
    834 
    835941  }
    836942
    … …  
    854960  G4double transitionEnergy = dedxCacheTransitionEnergy;
    855961
    856   if(iter != lossTableList.begin() && transitionEnergy < kineticEnergy) {
     962  if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
    857963     chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
    858964                                                                material,
    … …  
    862968     eloss *= chargeSquareRatioCorr;
    863969  }
    864   else if (iter == lossTableList.begin()) {
     970  else if (iter == lossTableList.end()) {
    865971
    866972     chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
    … …  
    875981  // overwrite the energy loss (i.e. when the effective charge approach is
    876982  // used)
    877   if(iter == lossTableList.begin()) {
     983  if(iter == lossTableList.end()) {
    878984
    879985     G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
    … …  
    9061012     particleChangeLoss -> ProposeNonIonizingEnergyDeposit(nloss);
    9071013  }
    908 
    909 }
    910 
     1014}
     1015
     1016// #########################################################################
    9111017
    9121018void G4IonParametrisedLossModel::BuildRangeVector(
    … …  
    10181124
    10191125  IonMatCouple ionMatCouple = std::make_pair(particle, material);
    1020  
     1126
    10211127  E[ionMatCouple] = energyRangeVector;
    10221128  r[ionMatCouple] = rangeEnergyVector;
    10231129}
    10241130
     1131// #########################################################################
     1132
     1133G4double G4IonParametrisedLossModel::GetRange(
     1134                  const G4ParticleDefinition* particle, // Projectile
     1135                  const G4Material* material,           // Target Material
     1136                  G4double kineticEnergy) {
     1137
     1138  G4double range = 0.0;
     1139
     1140  IonMatCouple couple = std::make_pair(particle, material);
     1141
     1142  EnergyRangeTable::iterator iter = E.find(couple);
     1143
     1144  if(iter == E.end()) {
     1145     G4cerr << "G4IonParametrisedLossModel::GetRange() No range vector found."
     1146            << G4endl;
     1147
     1148     G4cout << "   Ion-material pair: " << particle ->GetParticleName()
     1149            << "  " << material -> GetName()
     1150            << G4endl
     1151            << "   Available couples:"
     1152            << G4endl;
     1153
     1154     EnergyRangeTable::iterator iter_beg = E.begin();
     1155     EnergyRangeTable::iterator iter_end = E.end();
     1156
     1157     for(;iter_beg != iter_end; iter_beg++) {
     1158         IonMatCouple key = (*iter_beg).first;
     1159 
     1160         G4cout << "           " << (key.first) -> GetParticleName()
     1161                << "  " << (key.second) -> GetName()
     1162                << G4endl;
     1163     }
     1164  }
     1165  else {
     1166     G4PhysicsVector* energyRange = (*iter).second;
     1167
     1168     if(energyRange != 0) {
     1169        G4bool b;
     1170
     1171        // Computing range for kinetic energy:
     1172        range = energyRange -> GetValue(kineticEnergy, b);
     1173     }
     1174  }
     1175
     1176  return range;
     1177}
     1178
     1179// #########################################################################
    10251180
    10261181G4double G4IonParametrisedLossModel::ComputeLossForStep(
    … …  
    10611216        if(loss < 0.0) loss = 0.0;
    10621217     }
    1063 
    1064 #ifdef PRINT_DEBUG
    1065      G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() E = "
    1066             << kineticEnergy / MeV << " MeV * "
    1067             << value.energyScaling << " = "
    1068             << kineticEnergy * value.energyScaling / MeV
    1069             << " MeV, dE/dx = " << dedx / MeV * cm << " MeV/cm = "
    1070             << dedx/factor/MeV*cm << " * " << factor << " MeV/cm; index = "
    1071             << value.dEdxIndex << ", material = " << material -> GetName()
    1072             << G4endl;
    1073 #endif
    1074 
    10751218  }
    10761219
    10771220  return loss;
    10781221}
     1222
     1223// #########################################################################
     1224
     1225G4bool G4IonParametrisedLossModel::AddDEDXTable(
     1226                                const G4String& name,
     1227                                G4VIonDEDXTable* table,
     1228                                G4VIonDEDXScalingAlgorithm* algorithm) {
     1229
     1230  if(table == 0) {
     1231     G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
     1232            << " add table: Invalid pointer."
     1233            << G4endl;     
     1234
     1235     return false;
     1236  }
     1237
     1238  // Checking uniqueness of name
     1239  LossTableList::iterator iter = lossTableList.begin();
     1240  LossTableList::iterator iter_end = lossTableList.end();
     1241 
     1242  for(;iter != iter_end; iter++) {
     1243     G4String tableName = (*iter) -> GetName();
     1244
     1245     if(tableName == name) {
     1246        G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
     1247               << " add table: Name already exists."     
     1248               << G4endl;
     1249
     1250        return false;
     1251     }
     1252  }
     1253
     1254  G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
     1255  if(scalingAlgorithm == 0)
     1256     scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
     1257 
     1258  G4IonDEDXHandler* handler =
     1259                      new G4IonDEDXHandler(table, scalingAlgorithm, name);
     1260
     1261  lossTableList.push_front(handler);
     1262
     1263  return true;
     1264}
     1265
     1266// #########################################################################
     1267
     1268G4bool G4IonParametrisedLossModel::RemoveDEDXTable(
     1269                                 const G4String& name) {
     1270
     1271  LossTableList::iterator iter = lossTableList.begin();
     1272  LossTableList::iterator iter_end = lossTableList.end();
     1273 
     1274  for(;iter != iter_end; iter++) {
     1275     G4String tableName = (*iter) -> GetName();
     1276
     1277     if(tableName == name) {
     1278        delete (*iter);
     1279
     1280        lossTableList.erase(iter);
     1281        return true;
     1282     }
     1283  }
     1284
     1285  return false;
     1286}
     1287
     1288// #########################################################################
     1289

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermoreComptonModel.cc,v 1.1 2008/10/30 14:17:46 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
    28 //
     26// $Id: G4LivermoreComptonModel.cc,v 1.6 2009/04/18 18:29:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29//
     30// Author: Sebastien Inserti
     31//         30 October 2008
     32//
     33// History:
     34// --------
     35// 18 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     36//                  - apply internal high-energy limit only in constructor
     37//                  - do not apply low-energy limit (default is 0)
     38//                  - remove GetMeanFreePath method and table
     39//                  - added protection against numerical problem in energy sampling
     40//                  - use G4ElementSelector
    2941
    3042#include "G4LivermoreComptonModel.hh"
    … …  
    3749
    3850G4LivermoreComptonModel::G4LivermoreComptonModel(const G4ParticleDefinition*,
    39                                              const G4String& nam)
    40 :G4VEmModel(nam),isInitialised(false)
     51                                                 const G4String& nam)
     52  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
     53   scatterFunctionData(0),crossSectionHandler(0)
    4154{
    42   lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
     55  lowEnergyLimit = 250 * eV;
    4356  highEnergyLimit = 100 * GeV;
    44   SetLowEnergyLimit(lowEnergyLimit);
     57  //  SetLowEnergyLimit(lowEnergyLimit);
    4558  SetHighEnergyLimit(highEnergyLimit);
    4659
    47   verboseLevel= 0;
     60  verboseLevel=0 ;
    4861  // Verbosity scale:
    4962  // 0 = nothing
    … …  
    5265  // 3 = calculation of cross sections, file openings, sampling of atoms
    5366  // 4 = entering in methods
    54  
    55   G4cout << "Livermore Compton model is constructed " << G4endl
    56          << "Energy range: "
    57          << lowEnergyLimit / keV << " keV - "
    58          << highEnergyLimit / GeV << " GeV"
    59          << G4endl;
    60  
     67
     68  if(  verboseLevel>0 ) {
     69    G4cout << "Livermore Compton model is constructed " << G4endl
     70           << "Energy range: "
     71           << lowEnergyLimit / eV << " eV - "
     72           << highEnergyLimit / GeV << " GeV"
     73           << G4endl;
     74  }
    6175}
    6276
    … …  
    6579G4LivermoreComptonModel::~G4LivermoreComptonModel()
    6680{ 
    67   delete meanFreePathTable;
    68   delete crossSectionHandler;
    69   delete scatterFunctionData;
     81  if (crossSectionHandler) delete crossSectionHandler;
     82  if (scatterFunctionData) delete scatterFunctionData;
    7083}
    7184
    … …  
    7386
    7487void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
    75                                       const G4DataVector& cuts)
     88                                        const G4DataVector& cuts)
    7689{
    7790  if (verboseLevel > 3)
    7891    G4cout << "Calling G4LivermoreComptonModel::Initialise()" << G4endl;
    7992
    80   InitialiseElementSelectors(particle,cuts);
    81 
    82   // Energy limits
    83  
    84   if (LowEnergyLimit() < lowEnergyLimit)
    85     {
    86       G4cout << "G4LivermoreComptonModel: low energy limit increased from " <<
    87         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
    88       SetLowEnergyLimit(lowEnergyLimit);
    89     }
    90 
    91   if (HighEnergyLimit() > highEnergyLimit)
    92     {
    93       G4cout << "G4LivermoreComptonModel: high energy limit decreased from " <<
    94         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
    95       SetHighEnergyLimit(highEnergyLimit);
    96     }
    97 
     93  if (crossSectionHandler)
     94  {
     95    crossSectionHandler->Clear();
     96    delete crossSectionHandler;
     97  }
     98 
    9899  // Reading of data files - all materials are read
    99100 
    … …  
    113114  shellData.LoadData(file);
    114115
    115   meanFreePathTable = 0;
    116   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
    117  
    118116  if (verboseLevel > 2)
    119117    G4cout << "Loaded cross section files for Livermore Compton model" << G4endl;
    120118
    121   G4cout << "Livermore Compton model is initialized " << G4endl
    122          << "Energy range: "
    123          << LowEnergyLimit() / keV << " keV - "
    124          << HighEnergyLimit() / GeV << " GeV"
    125          << G4endl;
    126 
    127   //
    128  
     119  InitialiseElementSelectors(particle,cuts);
     120
     121  if(  verboseLevel>0 ) {
     122    G4cout << "Livermore Compton model is initialized " << G4endl
     123           << "Energy range: "
     124           << LowEnergyLimit() / eV << " eV - "
     125           << HighEnergyLimit() / GeV << " GeV"
     126           << G4endl;
     127  }
     128  // 
    129129  if(isInitialised) return;
    130 
    131   if(pParticleChange)
    132     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    133   else
    134     fParticleChange = new G4ParticleChangeForGamma();
    135    
     130  fParticleChange = GetParticleChangeForGamma();
    136131  isInitialised = true;
    137132}
    … …  
    148143    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << G4endl;
    149144
    150   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
     145  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
     146   
     147  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 
    151148  return cs;
    152149}
    … …  
    155152
    156153void G4LivermoreComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
    157                                               const G4MaterialCutsCouple* couple,
    158                                               const G4DynamicParticle* aDynamicGamma,
    159                                               G4double,
    160                                               G4double)
     154                                                const G4MaterialCutsCouple* couple,
     155                                                const G4DynamicParticle* aDynamicGamma,
     156                                                G4double, G4double)
    161157{
     158
    162159  // The scattered gamma energy is sampled according to Klein - Nishina formula.
    163160  // then accepted or rejected depending on the Scattering Function multiplied
    … …  
    174171  // (Nucl Phys 20(1960),15).
    175172
    176   if (verboseLevel > 3)
    177     G4cout << "Calling SampleSecondaries() of G4LivermoreComptonModel" << G4endl;
    178 
    179173  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
    180  
    181   if (photonEnergy0 <= lowEnergyLimit)
    182   {
     174
     175  if (verboseLevel > 3) {
     176    G4cout << "G4LivermoreComptonModel::SampleSecondaries() E(MeV)= "
     177           << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
     178           << G4endl;
     179  }
     180 
     181  // low-energy gamma is absorpted by this process
     182  if (photonEnergy0 <= lowEnergyLimit)
     183    {
    183184      fParticleChange->ProposeTrackStatus(fStopAndKill);
    184185      fParticleChange->SetProposedKineticEnergy(0.);
    185186      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
    186       // SI - IS THE FOLLOWING RETURN NECESSARY ?
    187187      return ;
    188   }
     188    }
    189189
    190190  G4double e0m = photonEnergy0 / electron_mass_c2 ;
    … …  
    192192
    193193  // Select randomly one element in the current material
    194   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     194  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     195  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     196  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
     197  G4int Z = (G4int)elm->GetZ();
    195198
    196199  G4double epsilon0 = 1. / (1. + 2. * e0m);
    … …  
    212215      if ( alpha1/(alpha1+alpha2) > G4UniformRand())
    213216      {
    214         epsilon = std::exp(-alpha1 * G4UniformRand());  // std::pow(epsilon0,G4UniformRand())
     217        // std::pow(epsilon0,G4UniformRand())
     218        epsilon = std::exp(-alpha1 * G4UniformRand()); 
    215219        epsilonSq = epsilon * epsilon;
    216220      }
    … …  
    238242  // Doppler broadening -  Method based on:
    239243  // Y. Namito, S. Ban and H. Hirayama,
    240   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
    241   // NIM A 349, pp. 489-494, 1994
     244  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
     245  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
    242246 
    243247  // Maximum number of sampling iterations
    … …  
    257261      eMax = photonEnergy0 - bindingE;
    258262     
    259       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
     263      // Randomly sample bound electron momentum
     264      // (memento: the data set is in Atomic Units)
    260265      G4double pSample = profileData.RandomSelectMomentum(Z,shell);
    261266      // Rescale from atomic units
    … …  
    299304
    300305  if (photonEnergy1 > 0.)
    301   {
    302     fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
    303   }
     306    {
     307      fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
     308    }
    304309  else
    305   {
    306     fParticleChange->SetProposedKineticEnergy(0.) ;
    307     fParticleChange->ProposeTrackStatus(fStopAndKill);
    308   }
     310    {
     311      photonEnergy1 = 0.;
     312      fParticleChange->SetProposedKineticEnergy(0.) ;
     313      fParticleChange->ProposeTrackStatus(fStopAndKill);   
     314    }
    309315
    310316  // Kinematics of the scattered electron
    311317  G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
     318
     319  // protection against negative final energy: no e- is created
     320  if(eKineticEnergy < 0.0) {
     321    fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0 - photonEnergy1);
     322    return;
     323  }
    312324  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
    313325
    … …  
    329341  eDirection.rotateUz(photonDirection0);
    330342
    331 // SI - The range test has been removed wrt original G4LowEnergyCompton class
     343  // SI - The range test has been removed wrt original G4LowEnergyCompton class
    332344
    333345  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
    334346 
    335   G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
     347  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
     348                                                 eDirection,eKineticEnergy) ;
    336349  fvect->push_back(dp);
    337350}
    338351
    339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    340 
    341 G4double G4LivermoreComptonModel::GetMeanFreePath(const G4Track& track,
    342                                              G4double, // previousStepSize
    343                                              G4ForceCondition*)
    344 {
    345   const G4DynamicParticle* photon = track.GetDynamicParticle();
    346   G4double energy = photon->GetKineticEnergy();
    347   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
    348   size_t materialIndex = couple->GetIndex();
    349 
    350   G4double meanFreePath;
    351   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
    352   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
    353   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
    354   return meanFreePath;
    355 }
    356 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreGammaConversionModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermoreGammaConversionModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
    28 //
     26// $Id: G4LivermoreGammaConversionModel.cc,v 1.6 2009/05/02 09:14:43 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29//
     30// Author: Sebastien Inserti
     31//         30 October 2008
     32//
     33// History:
     34// --------
     35// 12 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     36//                  - apply internal high-energy limit only in constructor
     37//                  - do not apply low-energy limit (default is 0)
     38//                  - use CLHEP electron mass for low-enegry limit
     39//                  - remove MeanFreePath method and table
     40
    2941
    3042#include "G4LivermoreGammaConversionModel.hh"
    … …  
    3749
    3850G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*,
    39                                              const G4String& nam)
    40 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false)
    41 {
    42   lowEnergyLimit = 1.022000 * MeV;
     51                                                                 const G4String& nam)
     52  :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),
     53   crossSectionHandler(0),meanFreePathTable(0)
     54{
     55  lowEnergyLimit = 2.0*electron_mass_c2;
    4356  highEnergyLimit = 100 * GeV;
    44  
    45   G4cout << "Livermore Gamma conversion is constructed " << G4endl
    46          << "Energy range: "
    47          << lowEnergyLimit / keV << " keV - "
    48          << highEnergyLimit / GeV << " GeV"
    49          << G4endl;
    50          
     57  SetHighEnergyLimit(highEnergyLimit);
     58         
    5159  verboseLevel= 0;
    5260  // Verbosity scale:
    … …  
    5765  // 4 = entering in methods
    5866
     67  if(verboseLevel > 0) {
     68    G4cout << "Livermore Gamma conversion is constructed " << G4endl
     69           << "Energy range: "
     70           << lowEnergyLimit / MeV << " MeV - "
     71           << highEnergyLimit / GeV << " GeV"
     72           << G4endl;
     73  }
    5974}
    6075
    … …  
    6378G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel()
    6479{ 
    65   delete meanFreePathTable;
    66   delete crossSectionHandler;
    67 }
    68 
    69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    70 
    71 void G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
    72                                        const G4DataVector& cuts)
     80  if (crossSectionHandler) delete crossSectionHandler;
     81}
     82
     83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     84
     85void
     86G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition*,
     87                                            const G4DataVector&)
    7388{
    7489  if (verboseLevel > 3)
    7590    G4cout << "Calling G4LivermoreGammaConversionModel::Initialise()" << G4endl;
    7691
    77   InitialiseElementSelectors(particle,cuts);
    78 
    79   // Energy limits
    80  
    81   if (LowEnergyLimit() < lowEnergyLimit)
     92  if (crossSectionHandler)
    8293  {
    83     G4cout << "G4LivermoreGammaConversionModel: low energy limit increased from " <<
    84         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
    85     SetLowEnergyLimit(lowEnergyLimit);
     94    crossSectionHandler->Clear();
     95    delete crossSectionHandler;
    8696  }
    8797
    88   if (HighEnergyLimit() > highEnergyLimit)
    89   {
    90     G4cout << "G4LivermoreGammaConversionModel: high energy limit decreased from " <<
    91         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
    92     SetHighEnergyLimit(highEnergyLimit);
    93   }
    94 
    9598  // Read data tables for all materials
    9699 
    97100  crossSectionHandler = new G4CrossSectionHandler();
    98   crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
     101  crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
    99102  G4String crossSectionFile = "pair/pp-cs-";
    100103  crossSectionHandler->LoadData(crossSectionFile);
    101104
    102   meanFreePathTable = 0;
    103   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
    104 
    105105  //
    106106 
    … …  
    108108    G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
    109109
    110   G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
    111          << "Energy range: "
    112          << LowEnergyLimit() / MeV << " MeV - "
    113          << HighEnergyLimit() / GeV << " GeV"
    114          << G4endl;
     110  if (verboseLevel > 0) {
     111    G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
     112           << "Energy range: "
     113           << LowEnergyLimit() / MeV << " MeV - "
     114           << HighEnergyLimit() / GeV << " GeV"
     115           << G4endl;
     116  }
    115117
    116118  if(isInitialised) return;
    117 
    118   if(pParticleChange)
    119     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    120   else
    121     fParticleChange = new G4ParticleChangeForGamma();
    122   isInitialised = true;}
    123 
    124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    125 
    126 G4double G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(
    127                                        const G4ParticleDefinition*,
    128                                              G4double GammaEnergy,
    129                                              G4double Z, G4double,
    130                                              G4double, G4double)
    131 {
    132   if (verboseLevel > 3)
    133     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel" << G4endl;
     119  fParticleChange = GetParticleChangeForGamma();
     120  isInitialised = true;
     121}
     122
     123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     124
     125G4double
     126G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     127                                                            G4double GammaEnergy,
     128                                                            G4double Z, G4double,
     129                                                            G4double, G4double)
     130{
     131  if (verboseLevel > 3) {
     132    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
     133           << G4endl;
     134  }
     135  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
    134136
    135137  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
    … …  
    173175    {
    174176      // Select randomly one element in the current material
    175       const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
     177      //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
     178      const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     179      const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
    176180
    177181      if (element == 0)
    178182        {
    179           G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0" << G4endl;
     183          G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
     184                 << G4endl;
     185          return;
    180186        }
    181187      G4IonisParamElm* ionisation = element->GetIonisation();
    182        if (ionisation == 0)
     188      if (ionisation == 0)
    183189        {
    184           G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0" << G4endl;
     190          G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
     191                 << G4endl;
     192          return;
    185193        }
    186194
    … …  
    272280  // distribution with respect to the Z axis along the parent photon
    273281 
    274 //  aParticleChange.SetNumberOfSecondaries(2) ;
    275282  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
    276283 
    277 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
     284  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
    278285
    279286  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
    … …  
    287294  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
    288295
    289 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
     296  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
    290297
    291298  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
    … …  
    337344}
    338345
    339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    340 
    341 G4double G4LivermoreGammaConversionModel::GetMeanFreePath(const G4Track& track,
    342                                                      G4double, // previousStepSize
    343                                                      G4ForceCondition*)
    344 {
    345   const G4DynamicParticle* photon = track.GetDynamicParticle();
    346   G4double energy = photon->GetKineticEnergy();
    347   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
    348   size_t materialIndex = couple->GetIndex();
    349 
    350   G4double meanFreePath;
    351   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
    352   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
    353   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
    354   return meanFreePath;
    355 }
    356 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermorePhotoElectricModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
    28 //
     26// $Id: G4LivermorePhotoElectricModel.cc,v 1.7 2009/04/18 18:29:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29//
     30// Author: Sebastien Inserti
     31//         30 October 2008
     32//
     33// History:
     34// --------
     35// 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     36//                  - apply internal high-energy limit only in constructor
     37//                  - do not apply low-energy limit (default is 0)
     38//                  - remove GetMeanFreePath method and table
     39//                  - simplify sampling of deexcitation by using cut in energy
     40//                  - added protection against numerical problem in energy sampling
     41//                  - use G4ElementSelector
    2942
    3043#include "G4LivermorePhotoElectricModel.hh"
    … …  
    3851G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4ParticleDefinition*,
    3952                                             const G4String& nam)
    40 :G4VEmModel(nam),isInitialised(false)
    41 {
    42   lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
     53:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
     54 crossSectionHandler(0),shellCrossSectionHandler(0),ElectronAngularGenerator(0)
     55{
     56  lowEnergyLimit = 250 * eV;
    4357  highEnergyLimit = 100 * GeV;
    44   SetLowEnergyLimit(lowEnergyLimit);
     58  //  SetLowEnergyLimit(lowEnergyLimit);
    4559  SetHighEnergyLimit(highEnergyLimit);
    46  
    47   G4cout << "Livermore Compton is constructed " << G4endl
    48          << "Energy range: "
    49          << lowEnergyLimit / keV << " keV - "
    50          << highEnergyLimit / GeV << " GeV"
    51          << G4endl;
    52  
     60
     61  ActivateAuger(false);
     62   
    5363  verboseLevel= 0;
    5464  // Verbosity scale:
    … …  
    5868  // 3 = calculation of cross sections, file openings, sampling of atoms
    5969  // 4 = entering in methods
     70  if(verboseLevel>0) {
     71    G4cout << "Livermore PhotoElectric is constructed " << G4endl
     72           << "Energy range: "
     73           << lowEnergyLimit / eV << " eV - "
     74           << highEnergyLimit / GeV << " GeV"
     75           << G4endl;
     76  }
    6077}
    6178
    … …  
    6481G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
    6582{ 
    66   delete meanFreePathTable;
    67   delete crossSectionHandler;
    68   delete shellCrossSectionHandler;
    69   delete ElectronAngularGenerator;
    70 }
    71 
    72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    73 
    74 void G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
    75                                        const G4DataVector& cuts)
     83  if (crossSectionHandler) delete crossSectionHandler;
     84  if (shellCrossSectionHandler) delete shellCrossSectionHandler;
     85  if (ElectronAngularGenerator) delete ElectronAngularGenerator;
     86}
     87
     88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     89
     90void
     91G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
     92                                          const G4DataVector&)
    7693{
    7794  if (verboseLevel > 3)
    7895    G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
    7996
    80   InitialiseElementSelectors(particle,cuts);
    81 
    82   // Energy limits
    83  
    84   if (LowEnergyLimit() < lowEnergyLimit)
     97  if (crossSectionHandler)
    8598  {
    86       G4cout << "G4LivermorePhotoElectricModel: low energy limit increased from " <<
    87         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" <<
    88         G4endl;
    89       SetLowEnergyLimit(lowEnergyLimit);
    90   }
    91  
    92   if (HighEnergyLimit() > highEnergyLimit)
     99    crossSectionHandler->Clear();
     100    delete crossSectionHandler;
     101  }
     102 
     103  if (shellCrossSectionHandler)
    93104  {
    94       G4cout << "G4LivermorePhotoElectricModel: high energy limit decreased from " <<
    95         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV"
    96              << G4endl;
    97       SetHighEnergyLimit(highEnergyLimit);
    98   }
    99 
     105    shellCrossSectionHandler->Clear();
     106    delete shellCrossSectionHandler;
     107  }
     108 
    100109  // Read data tables for all materials
    101110 
    … …  
    105114  crossSectionHandler->LoadData(crossSectionFile);
    106115
    107   meanFreePathTable = 0;
    108   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
    109 
    110116  shellCrossSectionHandler = new G4CrossSectionHandler();
    111117  shellCrossSectionHandler->Clear();
    … …  
    113119  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
    114120 
    115   // SI - Buggy default ?
     121  // SI - Simple generator is buggy
    116122  //generatorName = "geant4.6.2";
    117123  //ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator");              // default generator
    118 
    119   //
    120  
     124  ElectronAngularGenerator =
     125    new G4PhotoElectricAngularGeneratorSauterGavrila("GEANTSauterGavrilaGenerator");       
     126
     127  // 
    121128  if (verboseLevel > 2)
    122129    G4cout << "Loaded cross section files for Livermore PhotoElectric model" << G4endl;
    123130
    124   G4cout << "Livermore PhotoElectric model is initialized " << G4endl
    125          << "Energy range: "
    126          << LowEnergyLimit() / keV << " keV - "
    127          << HighEnergyLimit() / GeV << " GeV"
    128          << G4endl;
     131  //  InitialiseElementSelectors(particle,cuts);
     132
     133  if (verboseLevel > 0) {
     134    G4cout << "Livermore PhotoElectric model is initialized " << G4endl
     135           << "Energy range: "
     136           << LowEnergyLimit() / eV << " eV - "
     137           << HighEnergyLimit() / GeV << " GeV"
     138           << G4endl;
     139  }
    129140
    130141  if(isInitialised) return;
    131 
    132   if(pParticleChange)
    133     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    134   else
    135     fParticleChange = new G4ParticleChangeForGamma();
    136 
     142  fParticleChange = GetParticleChangeForGamma();
    137143  isInitialised = true;
    138144}
    … …  
    147153{
    148154  if (verboseLevel > 3)
    149     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel" << G4endl;
     155    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel"
     156           << G4endl;
     157
     158  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
     159    return 0;
    150160
    151161  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
    … …  
    155165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    156166
    157 void G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
    158                                               const G4MaterialCutsCouple* couple,
    159                                               const G4DynamicParticle* aDynamicGamma,
    160                                               G4double,
    161                                               G4double)
     167void
     168G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
     169                                                 const G4MaterialCutsCouple* couple,
     170                                                 const G4DynamicParticle* aDynamicGamma,
     171                                                 G4double,
     172                                                 G4double)
    162173{
    163174
    … …  
    172183  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
    173184 
     185  // kill incident photon
     186  fParticleChange->SetProposedKineticEnergy(0.);
     187  fParticleChange->ProposeTrackStatus(fStopAndKill);   
     188
     189  // low-energy gamma is absorpted by this process
    174190  if (photonEnergy <= lowEnergyLimit)
    175191    {
    176       fParticleChange->ProposeTrackStatus(fStopAndKill);
    177       fParticleChange->SetProposedKineticEnergy(0.);
    178192      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
    179       // SI - IS THE FOLLOWING RETURN NECESSARY ?
    180       return ;
     193      return;
    181194    }
    182195 
    183   G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection(); // Returns the normalized direction of the momentum
     196  // Returns the normalized direction of the momentum
     197  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
    184198
    185199  // Select randomly one element in the current material
    186   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
     200  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
     201  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     202  const G4Element* elm = SelectRandomAtom(couple->GetMaterial(),particle,photonEnergy);
     203  G4int Z = (G4int)elm->GetZ();
    187204
    188205  // Select the ionised shell in the current atom according to shell cross sections
    … …  
    195212  G4int shellId = shell->ShellId();
    196213
    197   // Create lists of pointers to DynamicParticles (photons and electrons)
    198   // (Is the electron vector necessary? To be checked)
    199   std::vector<G4DynamicParticle*>* photonVector = 0;
    200   std::vector<G4DynamicParticle*> electronVector;
    201 
    202   G4double energyDeposit = 0.0;
    203 
    204214  // Primary outcoming electron
    205215  G4double eKineticEnergy = photonEnergy - bindingEnergy;
    … …  
    209219  if (eKineticEnergy > 0.)
    210220    {
    211       // SI - Removed safety
    212      
    213       // Generate the electron only if with large enough range w.r.t. cuts and safety
    214       //G4double safety = aStep.GetPostStepPoint()->GetSafety();
    215 
    216       //if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
    217         {
    218 
    219           // Calculate direction of the photoelectron
    220           G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
    221           G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
    222 
    223           // The electron is created ...
    224           G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
    225                                                                electronDirection,
    226                                                                eKineticEnergy);
    227           electronVector.push_back(electron);
    228         }
    229       /*else
    230         {
    231           energyDeposit += eKineticEnergy;
    232         }*/
     221      // Calculate direction of the photoelectron
     222      G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
     223      G4ThreeVector electronDirection =
     224        ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,
     225                                                            eKineticEnergy,
     226                                                            gammaPolarization,
     227                                                            shellId);
     228
     229      // The electron is created ...
     230      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
     231                                                           electronDirection,
     232                                                           eKineticEnergy);
     233      fvect->push_back(electron);
    233234    }
    234235  else
    … …  
    237238    }
    238239
    239   G4int nElectrons = electronVector.size();
    240   size_t nTotPhotons = 0;
    241   G4int nPhotons=0;
    242   const G4ProductionCutsTable* theCoupleTable=
    243         G4ProductionCutsTable::GetProductionCutsTable();
    244 
    245   size_t index = couple->GetIndex();
    246   G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
    247   cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
    248  
    249   G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
    250   cute = std::min(cutForLowEnergySecondaryPhotons,cute);
    251 
    252   G4DynamicParticle* aPhoton;
    253 
    254   // Generation of fluorescence
    255   // Data in EADL are available only for Z > 5
    256   // Protection to avoid generating photons in the unphysical case of
    257   // shell binding energy > photon energy
    258   if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
    259     {
    260       photonVector = deexcitationManager.GenerateParticles(Z,shellId);
    261       nTotPhotons = photonVector->size();
    262       for (size_t k=0; k<nTotPhotons; k++)
    263         {
    264           aPhoton = (*photonVector)[k];
    265           if (aPhoton)
    266             {
    267               G4double itsCut = cutg;
    268               if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
    269               G4double itsEnergy = aPhoton->GetKineticEnergy();
    270 
    271               if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
    272                 {
    273                   nPhotons++;
    274                   // Local energy deposit is given as the sum of the
    275                   // energies of incident photons minus the energies
    276                   // of the outcoming fluorescence photons
    277                   bindingEnergy -= itsEnergy;
    278 
    279                 }
    280               else
    281                 {
    282                   delete aPhoton;
    283                   (*photonVector)[k] = 0;
    284                 }
    285             }
    286         }
    287     }
    288 
    289   energyDeposit += bindingEnergy;
    290 
    291   // Final state
    292  
    293   for (G4int l = 0; l<nElectrons; l++ )
    294     {
    295       aPhoton = electronVector[l];
    296       if(aPhoton) {
    297         fvect->push_back(aPhoton);
     240  // deexcitation
     241  if(DeexcitationFlag() && Z > 5) {
     242
     243    const G4ProductionCutsTable* theCoupleTable=
     244      G4ProductionCutsTable::GetProductionCutsTable();
     245
     246    size_t index = couple->GetIndex();
     247    G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
     248    G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
     249
     250    // Generation of fluorescence
     251    // Data in EADL are available only for Z > 5
     252    // Protection to avoid generating photons in the unphysical case of
     253    // shell binding energy > photon energy
     254    if (bindingEnergy > cutg || bindingEnergy > cute)
     255      {
     256        G4DynamicParticle* aPhoton;
     257        deexcitationManager.SetCutForSecondaryPhotons(cutg);
     258        deexcitationManager.SetCutForAugerElectrons(cute);
     259 
     260        std::vector<G4DynamicParticle*>* photonVector =
     261          deexcitationManager.GenerateParticles(Z,shellId);
     262        size_t nTotPhotons = photonVector->size();
     263        for (size_t k=0; k<nTotPhotons; k++)
     264          {
     265            aPhoton = (*photonVector)[k];
     266            if (aPhoton)
     267              {
     268                G4double itsEnergy = aPhoton->GetKineticEnergy();
     269                if (itsEnergy <= bindingEnergy)
     270                  {
     271                    // Local energy deposit is given as the sum of the
     272                    // energies of incident photons minus the energies
     273                    // of the outcoming fluorescence photons
     274                    bindingEnergy -= itsEnergy;
     275                    fvect->push_back(aPhoton);
     276                  }
     277                else
     278                  {
     279                    // abnormal case of energy non-conservation
     280                    delete aPhoton;
     281                    (*photonVector)[k] = 0;
     282                  }
     283              }
     284          }
     285        delete photonVector;
    298286      }
    299     }
    300   for ( size_t ll = 0; ll < nTotPhotons; ll++)
    301     {
    302       aPhoton = (*photonVector)[ll];
    303       if(aPhoton) {
    304         fvect->push_back(aPhoton);
    305       }
    306     }
    307 
    308   delete photonVector;
    309 
    310   if (energyDeposit < 0)
    311     {
    312       G4cout << "WARNING - "
    313              << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
    314              << G4endl;
    315       energyDeposit = 0;
    316     }
    317 
    318   // kill incident photon
    319   fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
    320   fParticleChange->SetProposedKineticEnergy(0.);
    321   fParticleChange->ProposeTrackStatus(fStopAndKill);   
    322   fParticleChange->ProposeLocalEnergyDeposit(energyDeposit);
     287  }
     288  // excitation energy left
     289  fParticleChange->ProposeLocalEnergyDeposit(bindingEnergy);
    323290}
    324291
    … …  
    348315//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    349316
    350 void G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* distribution)
    351 {
    352   ElectronAngularGenerator = distribution;
     317void
     318G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* dist)
     319{
     320  ElectronAngularGenerator = dist;
    353321  ElectronAngularGenerator->PrintGeneratorInformation();
    354322}
    … …  
    361329    {
    362330      delete ElectronAngularGenerator;
    363       ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
     331      ElectronAngularGenerator =
     332        new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
    364333      generatorName = name;
    365334    }
    … …  
    367336    {
    368337      delete ElectronAngularGenerator;
    369       ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
     338      ElectronAngularGenerator =
     339        new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
    370340      generatorName = name;
    371341    }
    … …  
    373343    {
    374344      delete ElectronAngularGenerator;
    375       ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
     345      ElectronAngularGenerator =
     346        new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
    376347      generatorName = name;
    377348    }
    … …  
    384355}
    385356
    386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    387 
    388 G4double G4LivermorePhotoElectricModel::GetMeanFreePath(const G4Track& track,
    389                                                    G4double, // previousStepSize
    390                                                G4ForceCondition*)
    391 {
    392   const G4DynamicParticle* photon = track.GetDynamicParticle();
    393   G4double energy = photon->GetKineticEnergy();
    394   G4Material* material = track.GetMaterial();
    395   //  size_t materialIndex = material->GetIndex();
    396 
    397   G4double meanFreePath = DBL_MAX;
    398 
    399   //  if (energy > highEnergyLimit)
    400   //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
    401   //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
    402   //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
    403 
    404   G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
    405   if(cross > 0.0) meanFreePath = 1.0/cross;
    406 
    407   return meanFreePath;
    408 }
    409 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermorePolarizedComptonModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LivermorePolarizedComptonModel.cc,v 1.6 2009/05/03 08:29:55 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
     29// History:
     30// --------
     31// 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
     32//
     33// Cleanup initialisation and generation of secondaries:
     34//                  - apply internal high-energy limit only in constructor
     35//                  - do not apply low-energy limit (default is 0)
     36//                  - remove GetMeanFreePath method and table
     37//                  - added protection against numerical problem in energy sampling
     38//                  - use G4ElementSelector
    2939
    3040#include "G4LivermorePolarizedComptonModel.hh"
    … …  
    3848G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
    3949                                             const G4String& nam)
    40 :G4VEmModel(nam),isInitialised(false)
    41 {
    42   lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
     50:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
     51{
     52  lowEnergyLimit = 250 * eV;
    4353  highEnergyLimit = 100 * GeV;
    44   SetLowEnergyLimit(lowEnergyLimit);
     54  //SetLowEnergyLimit(lowEnergyLimit);
    4555  SetHighEnergyLimit(highEnergyLimit);
    4656 
    … …  
    5363  // 4 = entering in methods
    5464
     65  if( verboseLevel>0 ) {
    5566  G4cout << "Livermore Polarized Compton is constructed " << G4endl
    5667         << "Energy range: "
    57          << lowEnergyLimit / keV << " keV - "
     68         << lowEnergyLimit / eV << " eV - "
    5869         << highEnergyLimit / GeV << " GeV"
    5970         << G4endl;
    60 
     71  }
    6172}
    6273
    … …  
    6576G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
    6677{ 
    67   delete meanFreePathTable;
    68   delete crossSectionHandler;
    69   delete scatterFunctionData;
     78  if (meanFreePathTable)   delete meanFreePathTable;
     79  if (crossSectionHandler) delete crossSectionHandler;
     80  if (scatterFunctionData) delete scatterFunctionData;
    7081}
    7182
    … …  
    7889    G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
    7990
    80   InitialiseElementSelectors(particle,cuts);
    81 
    82   // Energy limits
    83  
    84   if (LowEnergyLimit() < lowEnergyLimit)
    85     {
    86       G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " <<
    87         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
    88       SetLowEnergyLimit(lowEnergyLimit);
    89     }
    90 
    91   if (HighEnergyLimit() > highEnergyLimit)
    92     {
    93       G4cout << "G4LivermorePolarizedComptonModel: high energy limit decreased from " <<
    94         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
    95       SetHighEnergyLimit(highEnergyLimit);
    96     }
     91  if (crossSectionHandler)
     92  {
     93    crossSectionHandler->Clear();
     94    delete crossSectionHandler;
     95  }
    9796
    9897  // Reading of data files - all materials are read
    … …  
    116115  shellData.LoadData(file);
    117116
    118   //
    119117  if (verboseLevel > 2)
    120118    G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
    121119
    122   G4cout << "Livermore Polarized Compton model is initialized " << G4endl
     120  InitialiseElementSelectors(particle,cuts);
     121
     122  if(  verboseLevel>0 ) {
     123    G4cout << "Livermore Polarized Compton model is initialized " << G4endl
    123124         << "Energy range: "
    124          << LowEnergyLimit() / keV << " keV - "
     125         << LowEnergyLimit() / eV << " eV - "
    125126         << HighEnergyLimit() / GeV << " GeV"
    126127         << G4endl;
    127 
     128  }
     129 
    128130  //
    129131   
    130132  if(isInitialised) return;
    131 
    132   if(pParticleChange)
    133     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    134   else
    135     fParticleChange = new G4ParticleChangeForGamma();
    136    
     133  fParticleChange = GetParticleChangeForGamma();
    137134  isInitialised = true;
    138135}
    … …  
    149146    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
    150147
     148  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
     149
    151150  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
    152151  return cs;
    … …  
    203202      fParticleChange->SetProposedKineticEnergy(0.);
    204203      fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
    205       // SI - IS THE FOLLOWING RETURN NECESSARY ?
    206204      return;
    207205    }
    … …  
    210208
    211209  // Select randomly one element in the current material
    212 
    213   G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
     210  //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
     211  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     212  const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
     213  G4int Z = (G4int)elm->GetZ();
    214214
    215215  // Sample the energy and the polarization of the scattered photon
    … …  
    429429  else
    430430    {
     431      gammaEnergy1 = 0.;
    431432      fParticleChange->SetProposedKineticEnergy(0.) ;
    432433      fParticleChange->ProposeTrackStatus(fStopAndKill);
    … …  
    439440  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
    440441
     442  // SI -protection against negative final energy: no e- is created
     443  // like in G4LivermoreComptonModel.cc
     444  if(ElecKineEnergy < 0.0) {
     445    fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
     446    return;
     447  }
     448 
    441449  // SI - Removed range test
    442450 
    … …  
    664672}
    665673
    666 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    667 
    668 G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track,
    669                                                       G4double,
    670                                                       G4ForceCondition*)
    671 {
    672   const G4DynamicParticle* photon = track.GetDynamicParticle();
    673   G4double energy = photon->GetKineticEnergy();
    674   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
    675   size_t materialIndex = couple->GetIndex();
    676   G4double meanFreePath;
    677   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
    678   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
    679   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
    680   return meanFreePath;
    681 }
    682 
    683 
     674

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedRayleighModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermorePolarizedRayleighModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LivermorePolarizedRayleighModel.cc,v 1.5 2009/05/02 15:20:53 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
     29// History:
     30// --------
     31// 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreRayleighModel.cc
     32//
     33// Cleanup initialisation and generation of secondaries:
     34//                  - apply internal high-energy limit only in constructor
     35//                  - do not apply low-energy limit (default is 0)
     36//                  - remove GetMeanFreePath method and table
     37//                  - remove initialisation of element selector
     38//                  - use G4ElementSelector
    2939
    3040#include "G4LivermorePolarizedRayleighModel.hh"
    … …  
    3848G4LivermorePolarizedRayleighModel::G4LivermorePolarizedRayleighModel(const G4ParticleDefinition*,
    3949                                             const G4String& nam)
    40 :G4VEmModel(nam),isInitialised(false)
    41 {
    42   lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
     50:G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),formFactorData(0)
     51{
     52  lowEnergyLimit = 250 * eV;
    4353  highEnergyLimit = 100 * GeV;
    4454 
    45   SetLowEnergyLimit(lowEnergyLimit);
     55  //SetLowEnergyLimit(lowEnergyLimit);
    4656  SetHighEnergyLimit(highEnergyLimit);
    4757  //
    … …  
    5464  // 4 = entering in methods
    5565
    56   G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
     66  if(verboseLevel > 0) {
     67    G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
    5768         << "Energy range: "
    58          << lowEnergyLimit / keV << " keV - "
     69         << lowEnergyLimit / eV << " eV - "
    5970         << highEnergyLimit / GeV << " GeV"
    6071         << G4endl;
     72  }
    6173}
    6274
    … …  
    6577G4LivermorePolarizedRayleighModel::~G4LivermorePolarizedRayleighModel()
    6678{ 
    67   delete crossSectionHandler;
    68   delete formFactorData;
     79  if (crossSectionHandler) delete crossSectionHandler;
     80  if (formFactorData) delete formFactorData;
    6981}
    7082
    … …  
    8496    G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
    8597
    86   InitialiseElementSelectors(particle,cuts);
    87 
    88   // Energy limits
    89  
    90   if (LowEnergyLimit() < lowEnergyLimit)
    91     {
    92       G4cout << "G4LivermorePolarizedRayleighModel: low energy limit increased from " <<
    93         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
    94       SetLowEnergyLimit(lowEnergyLimit);
    95     }
    96   if (HighEnergyLimit() > highEnergyLimit)
    97     {
    98       G4cout << "G4LivermorePolarizedRayleighModel: high energy limit decreased from " <<
    99         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
    100       SetHighEnergyLimit(highEnergyLimit);
    101     }
    102    
     98  if (crossSectionHandler)
     99  {
     100    crossSectionHandler->Clear();
     101    delete crossSectionHandler;
     102  }
     103 
    103104  // Read data files for all materials
    104105
    … …  
    113114  formFactorData->LoadData(formFactorFile);
    114115
     116  InitialiseElementSelectors(particle,cuts);
     117
    115118  //
    116119  if (verboseLevel > 2)
    117120    G4cout << "Loaded cross section files for Livermore Polarized Rayleigh model" << G4endl;
    118121
    119   G4cout << "Livermore Polarized Rayleigh model is initialized " << G4endl
     122  InitialiseElementSelectors(particle,cuts);
     123
     124  if (verboseLevel > 0) {
     125    G4cout << "Livermore Polarized Rayleigh model is initialized " << G4endl
    120126         << "Energy range: "
    121          << LowEnergyLimit() / keV << " keV - "
     127         << LowEnergyLimit() / eV << " eV - "
    122128         << HighEnergyLimit() / GeV << " GeV"
    123129         << G4endl;
     130         }
    124131
    125132  if(isInitialised) return;
    126 
    127   if(pParticleChange)
    128     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    129   else
    130     fParticleChange = new G4ParticleChangeForGamma();
    131    
     133  fParticleChange = GetParticleChangeForGamma();
    132134  isInitialised = true;
    133135}
    … …  
    144146    G4cout << "Calling CrossSectionPerAtom() of G4LivermorePolarizedRayleighModel" << G4endl;
    145147
     148  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
     149
    146150  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
    147151  return cs;
    … …  
    166170      fParticleChange->SetProposedKineticEnergy(0.);
    167171      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
    168       // SI - IS THE FOLLOWING RETURN NECESSARY ?
    169172      return ;
    170173  }
    … …  
    173176
    174177  // Select randomly one element in the current material
    175   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     178  // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     179  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     180  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
     181  G4int Z = (G4int)elm->GetZ();
    176182
    177183  G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreRayleighModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermoreRayleighModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
    28 //
     26// $Id: G4LivermoreRayleighModel.cc,v 1.6 2009/04/18 18:29:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     28//
     29// Author: Sebastien Inserti
     30//         30 October 2008
     31//
     32// History:
     33// --------
     34// 18 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     35//                  - apply internal high-energy limit only in constructor
     36//                  - do not apply low-energy limit (default is 0)
     37//                  - remove GetMeanFreePath method and table
     38//                  - remove initialisation of element selector
     39//                  - use G4ElementSelector
    2940
    3041#include "G4LivermoreRayleighModel.hh"
    … …  
    3748
    3849G4LivermoreRayleighModel::G4LivermoreRayleighModel(const G4ParticleDefinition*,
    39                                              const G4String& nam)
    40 :G4VEmModel(nam),isInitialised(false)
    41 {
    42   lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
     50                                                   const G4String& nam)
     51  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
     52   formFactorData(0),crossSectionHandler(0)
     53{
     54  lowEnergyLimit = 250 * eV;
    4355  highEnergyLimit = 100 * GeV;
    4456 
    45   SetLowEnergyLimit(lowEnergyLimit);
     57  //  SetLowEnergyLimit(lowEnergyLimit);
    4658  SetHighEnergyLimit(highEnergyLimit);
    4759  //
    … …  
    5466  // 4 = entering in methods
    5567
    56   G4cout << "Livermore Rayleigh is constructed " << G4endl
    57          << "Energy range: "
    58          << lowEnergyLimit / keV << " keV - "
    59          << highEnergyLimit / GeV << " GeV"
    60          << G4endl;
     68  if(verboseLevel > 0) {
     69    G4cout << "Livermore Rayleigh is constructed " << G4endl
     70           << "Energy range: "
     71           << lowEnergyLimit / eV << " eV - "
     72           << highEnergyLimit / GeV << " GeV"
     73           << G4endl;
     74  }
    6175}
    6276
    … …  
    6579G4LivermoreRayleighModel::~G4LivermoreRayleighModel()
    6680{ 
    67   delete meanFreePathTable;
    68   delete crossSectionHandler;
    69   delete formFactorData;
     81  if (crossSectionHandler) delete crossSectionHandler;
     82  if (formFactorData) delete formFactorData;
    7083}
    7184
    … …  
    7386
    7487void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle,
    75                                        const G4DataVector& cuts)
     88                                          const G4DataVector& cuts)
    7689{
    7790  if (verboseLevel > 3)
    7891    G4cout << "Calling G4LivermoreRayleighModel::Initialise()" << G4endl;
    7992
    80   InitialiseElementSelectors(particle,cuts);
    81 
    82   // Energy limits
     93  if (crossSectionHandler)
     94  {
     95    crossSectionHandler->Clear();
     96    delete crossSectionHandler;
     97  }
    8398 
    84   if (LowEnergyLimit() < lowEnergyLimit)
    85   {
    86     G4cout << "G4LivermoreRayleighModel: low energy limit increased from " <<
    87         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
    88     SetLowEnergyLimit(lowEnergyLimit);
    89   }
    90 
    91   if (HighEnergyLimit() > highEnergyLimit)
    92   {
    93     G4cout << "G4LivermoreRayleighModel: high energy limit decreased from " <<
    94         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
    95     SetHighEnergyLimit(highEnergyLimit);
    96   }
    97 
    9899  // Data are read for all materials
    99100 
    … …  
    103104  crossSectionHandler->LoadData(crossSectionFile);
    104105
    105   meanFreePathTable = 0;
    106   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
    107 
    108106  G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
    109107  G4String formFactorFile = "rayl/re-ff-";
    … …  
    111109  formFactorData->LoadData(formFactorFile);
    112110
    113   //
    114  
     111  InitialiseElementSelectors(particle,cuts);
     112
     113  // 
    115114  if (verboseLevel > 2)
    116115    G4cout << "Loaded cross section files for Livermore Rayleigh model" << G4endl;
    117116
    118   G4cout << "Livermore Rayleigh model is initialized " << G4endl
    119          << "Energy range: "
    120          << LowEnergyLimit() / keV << " keV - "
    121          << HighEnergyLimit() / GeV << " GeV"
    122          << G4endl;
     117  if (verboseLevel > 0) {
     118    G4cout << "Livermore Rayleigh model is initialized " << G4endl
     119           << "Energy range: "
     120           << LowEnergyLimit() / eV << " eV - "
     121           << HighEnergyLimit() / GeV << " GeV"
     122           << G4endl;
     123  }
    123124
    124125  if(isInitialised) return;
    125 
    126   if(pParticleChange)
    127     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    128   else
    129     fParticleChange = new G4ParticleChangeForGamma();
    130 
     126  fParticleChange = GetParticleChangeForGamma();
    131127  isInitialised = true;
    132128
    … …  
    144140    G4cout << "Calling CrossSectionPerAtom() of G4LivermoreRayleighModel" << G4endl;
    145141
     142  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
     143    return 0.0;
     144
    146145  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
    147146  return cs;
    … …  
    160159
    161160  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
    162  
     161
     162  // absorption of low-energy gamma 
    163163  if (photonEnergy0 <= lowEnergyLimit)
    164   {
     164    {
    165165      fParticleChange->ProposeTrackStatus(fStopAndKill);
    166166      fParticleChange->SetProposedKineticEnergy(0.);
    167167      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
    168       // SI - IS THE FOLLOWING RETURN NECESSARY ?
    169168      return ;
    170   }
     169    }
    171170
    172171  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
    173172
    174173  // Select randomly one element in the current material
    175   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     174  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
     175  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
     176  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
     177  G4int Z = (G4int)elm->GetZ();
    176178
    177179  // Sample the angle of the scattered photon
    … …  
    188190    {
    189191      do
    190       {
    191       cosTheta = 2. * G4UniformRand() - 1.;
    192       fcostheta = ( 1. + cosTheta*cosTheta)/2.;
    193       } while (fcostheta < G4UniformRand());
     192        {
     193          cosTheta = 2. * G4UniformRand() - 1.;
     194          fcostheta = ( 1. + cosTheta*cosTheta)/2.;
     195        } while (fcostheta < G4UniformRand());
    194196
    195197      G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
    196198      x = sinThetaHalf / (wlPhoton/cm);
    197199      if (x > 1.e+005)
    198          dataFormFactor = formFactorData->FindValue(x,Z-1);
     200        {
     201          dataFormFactor = formFactorData->FindValue(x,Z-1);
     202        }
    199203      else
    200          dataFormFactor = formFactorData->FindValue(0.,Z-1);
     204        {
     205          dataFormFactor = formFactorData->FindValue(0.,Z-1);
     206        }
    201207      randomFormFactor = G4UniformRand() * Z * Z;
    202208      sinTheta = std::sqrt(1. - cosTheta*cosTheta);
    … …  
    219225}
    220226
    221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    222 
    223 G4double G4LivermoreRayleighModel::GetMeanFreePath(const G4Track& track,
    224                                               G4double, // previousStepSize
    225                                               G4ForceCondition*)
    226 {
    227   const G4DynamicParticle* photon = track.GetDynamicParticle();
    228   G4double energy = photon->GetKineticEnergy();
    229   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
    230   size_t materialIndex = couple->GetIndex();
    231 
    232   G4double meanFreePath;
    233   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
    234   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
    235   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
    236   return meanFreePath;
    237 }
    238 
     227

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyBremsstrahlung.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LowEnergyBremsstrahlung.cc,v 1.71 2006/06/29 19:40:13 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LowEnergyBremsstrahlung.cc,v 1.72 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    9393//  angularDistribution->PrintGeneratorInformation();
    9494  TsaiAngularDistribution = new G4ModifiedTsai("TsaiGenerator");
     95
     96
     97   G4cout << G4endl;
     98   G4cout << "*******************************************************************************" << G4endl;
     99   G4cout << "*******************************************************************************" << G4endl;
     100   G4cout << "   The class G4LowEnergyBremsstrahlung is NOT SUPPORTED ANYMORE. " << G4endl;
     101   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     102   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     103   G4cout << "*******************************************************************************" << G4endl;
     104   G4cout << "*******************************************************************************" << G4endl;
     105   G4cout << G4endl;
    95106}
    96107

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyCompton.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LowEnergyCompton.cc,v 1.47 2008/12/18 13:01:28 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LowEnergyCompton.cc,v 1.48 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: A. Forti
    … …  
    104104              << G4endl;
    105105     }
     106
     107   G4cout << G4endl;
     108   G4cout << "*******************************************************************************" << G4endl;
     109   G4cout << "*******************************************************************************" << G4endl;
     110   G4cout << "   The class G4LowEnergyCompton is NOT SUPPORTED ANYMORE. " << G4endl;
     111   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     112   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     113   G4cout << "*******************************************************************************" << G4endl;
     114   G4cout << "*******************************************************************************" << G4endl;
     115   G4cout << G4endl;
    106116}
    107117

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyGammaConversion.cc

    r1007 r1055  
    2626// --------------------------------------------------------------------
    2727///
    28 // $Id: G4LowEnergyGammaConversion.cc,v 1.36 2006/06/29 19:40:17 gunter Exp $
    29 // GEANT4 tag $Name: geant4-09-02 $
     28// $Id: G4LowEnergyGammaConversion.cc,v 1.37 2009/05/02 09:59:16 sincerti Exp $
     29// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    3030//
    3131//
    … …  
    104104              << G4endl;
    105105     }
     106
     107   G4cout << G4endl;
     108   G4cout << "*******************************************************************************" << G4endl;
     109   G4cout << "*******************************************************************************" << G4endl;
     110   G4cout << "   The class G4LowEnergyGammaConversion is NOT SUPPORTED ANYMORE. " << G4endl;
     111   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     112   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     113   G4cout << "*******************************************************************************" << G4endl;
     114   G4cout << "*******************************************************************************" << G4endl;
     115   G4cout << G4endl;
    106116}
    107117 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyIonisation.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LowEnergyIonisation.cc,v 1.103 2008/05/02 19:23:38 pia Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LowEnergyIonisation.cc,v 1.104 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    131131  cutForElectrons = 250.0*eV;
    132132  verboseLevel = 0;
     133
     134   G4cout << G4endl;
     135   G4cout << "*******************************************************************************" << G4endl;
     136   G4cout << "*******************************************************************************" << G4endl;
     137   G4cout << "   The class G4LowEnergyIonisation is NOT SUPPORTED ANYMORE. " << G4endl;
     138   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     139   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     140   G4cout << "*******************************************************************************" << G4endl;
     141   G4cout << "*******************************************************************************" << G4endl;
     142   G4cout << G4endl;
    133143}
    134144

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyPhotoElectric.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4LowEnergyPhotoElectric.cc,v 1.56 2006/06/29 19:40:23 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4LowEnergyPhotoElectric.cc,v 1.57 2009/05/02 09:59:16 sincerti Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Author: A. Forti
    … …  
    122122             << G4endl;
    123123    }
     124
     125   G4cout << G4endl;
     126   G4cout << "*******************************************************************************" << G4endl;
     127   G4cout << "*******************************************************************************" << G4endl;
     128   G4cout << "   The class G4LowEnergyPhotoElectric is NOT SUPPORTED ANYMORE. " << G4endl;
     129   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     130   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     131   G4cout << "*******************************************************************************" << G4endl;
     132   G4cout << "*******************************************************************************" << G4endl;
     133   G4cout << G4endl;
    124134}
    125135

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyPolarizedCompton.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4LowEnergyPolarizedCompton.cc,v 1.25 2008/05/02 19:23:38 pia Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4LowEnergyPolarizedCompton.cc,v 1.26 2009/05/02 09:59:16 sincerti Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// ------------------------------------------------------------
    … …  
    122122              << G4endl;
    123123     }
     124
     125   G4cout << G4endl;
     126   G4cout << "*******************************************************************************" << G4endl;
     127   G4cout << "*******************************************************************************" << G4endl;
     128   G4cout << "   The class G4LowEnergyPolarizedCompton is NOT SUPPORTED ANYMORE. " << G4endl;
     129   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     130   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     131   G4cout << "*******************************************************************************" << G4endl;
     132   G4cout << "*******************************************************************************" << G4endl;
     133   G4cout << G4endl;
    124134}
    125135

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyPolarizedRayleigh.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LowEnergyPolarizedRayleigh.cc,v 1.7 2006/06/29 19:40:27 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LowEnergyPolarizedRayleigh.cc,v 1.8 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    6363 
    6464 
     65   G4cout << G4endl;
     66   G4cout << "*******************************************************************************" << G4endl;
     67   G4cout << "*******************************************************************************" << G4endl;
     68   G4cout << "   The class G4LowEnergyPolarizedRayleigh is NOT SUPPORTED ANYMORE. " << G4endl;
     69   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     70   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     71   G4cout << "*******************************************************************************" << G4endl;
     72   G4cout << "*******************************************************************************" << G4endl;
     73   G4cout << G4endl;
    6574}
    6675

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyRayleigh.cc

    r1007 r1055  
    2626// --------------------------------------------------------------------
    2727//
    28 // $Id: G4LowEnergyRayleigh.cc,v 1.37 2006/06/29 19:40:29 gunter Exp $
    29 // GEANT4 tag $Name: geant4-09-02 $
     28// $Id: G4LowEnergyRayleigh.cc,v 1.38 2009/05/02 09:59:16 sincerti Exp $
     29// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    3030//
    3131// Author: A. Forti
    … …  
    9999              << G4endl;
    100100     }
     101
     102   G4cout << G4endl;
     103   G4cout << "*******************************************************************************" << G4endl;
     104   G4cout << "*******************************************************************************" << G4endl;
     105   G4cout << "   The class G4LowEnergyRayleigh is NOT SUPPORTED ANYMORE. " << G4endl;
     106   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     107   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     108   G4cout << "*******************************************************************************" << G4endl;
     109   G4cout << "*******************************************************************************" << G4endl;
     110   G4cout << G4endl;
    101111}
    102112

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilation.cc

    r819 r1055  
    6464             << G4endl;
    6565    }
     66
     67   G4cout << G4endl;
     68   G4cout << "*******************************************************************************" << G4endl;
     69   G4cout << "*******************************************************************************" << G4endl;
     70   G4cout << "   The class G4PenelopeAnnihilation is NOT SUPPORTED ANYMORE. " << G4endl;
     71   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     72   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     73   G4cout << "*******************************************************************************" << G4endl;
     74   G4cout << "*******************************************************************************" << G4endl;
     75   G4cout << G4endl;
    6676}
    6777

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeAnnihilationModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeAnnihilationModel.cc,v 1.3 2009/04/17 10:29:20 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    3232// --------
    3333// 29 Oct 2008   L Pandola    Migration from process to model
    34 //
     34// 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     35//                  - apply internal high-energy limit only in constructor
     36//                  - do not apply low-energy limit (default is 0)
     37//                  - do not use G4ElementSelector
    3538
    3639#include "G4PenelopeAnnihilationModel.hh"
    … …  
    4851  :G4VEmModel(nam),isInitialised(false)
    4952{
    50   fIntrinsicLowEnergyLimit = 0.0*eV;
     53  fIntrinsicLowEnergyLimit = 0.0;
    5154  fIntrinsicHighEnergyLimit = 100.0*GeV;
    52   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     55  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    5356  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    5457 
    … …  
    7376//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    7477
    75 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition* particle,
    76                                        const G4DataVector& cuts)
     78void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*,
     79                                             const G4DataVector&)
    7780{
    7881  if (verboseLevel > 3)
    7982    G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
    8083
    81   InitialiseElementSelectors(particle,cuts);
    82   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    83     {
    84       G4cout << "G4PenelopeAnnihilationModel: low energy limit increased from " <<
    85         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
    86       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    87     }
    88   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    89     {
    90       G4cout << "G4PenelopeAnnihilationModel: high energy limit decreased from " <<
    91         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
    92       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    93     }
    94 
    95   G4cout << "Penelope Annihilation model is initialized " << G4endl
    96          << "Energy range: "
    97          << LowEnergyLimit() / keV << " keV - "
    98          << HighEnergyLimit() / GeV << " GeV"
    99          << G4endl;
     84  if(verboseLevel > 0) {
     85    G4cout << "Penelope Annihilation model is initialized " << G4endl
     86           << "Energy range: "
     87           << LowEnergyLimit() / keV << " keV - "
     88           << HighEnergyLimit() / GeV << " GeV"
     89           << G4endl;
     90  }
    10091
    10192  if(isInitialised) return;
    102 
    103   if(pParticleChange)
    104     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    105   else
    106     fParticleChange = new G4ParticleChangeForGamma();
     93  fParticleChange = GetParticleChangeForGamma();
    10794  isInitialised = true;
    10895}
    … …  
    155142
    156143  G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
    157  
    158   if (kineticEnergy == 0)
     144
     145  // kill primary
     146  fParticleChange->SetProposedKineticEnergy(0.);
     147  fParticleChange->ProposeTrackStatus(fStopAndKill);
     148 
     149  if (kineticEnergy == 0.0)
    159150    {
    160151      //Old AtRestDoIt
    … …  
    170161      fvect->push_back(firstGamma);
    171162      fvect->push_back(secondGamma);
    172       fParticleChange->SetProposedKineticEnergy(0.);
    173       fParticleChange->ProposeTrackStatus(fStopAndKill);
    174163      return;
    175164    }
    … …  
    227216                                                           photon2Energy);
    228217  fvect->push_back(aParticle2);
    229   fParticleChange->SetProposedKineticEnergy(0.);
    230   fParticleChange->ProposeTrackStatus(fStopAndKill);
    231218
    232219  if (verboseLevel > 1)
    … …  
    244231    }
    245232  if (verboseLevel > 0)
    246     {
    247      
     233    {     
    248234      G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
    249235      if (energyDiff > 0.05*keV)

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlung.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeBremsstrahlung.cc,v 1.18 2006/06/29 19:40:35 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeBremsstrahlung.cc,v 1.19 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    7171  cutForPhotons = 0.;
    7272  verboseLevel = 0;
     73
     74   G4cout << G4endl;
     75   G4cout << "*******************************************************************************" << G4endl;
     76   G4cout << "*******************************************************************************" << G4endl;
     77   G4cout << "   The class G4PenelopeBremsstrahlung is NOT SUPPORTED ANYMORE. " << G4endl;
     78   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     79   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     80   G4cout << "*******************************************************************************" << G4endl;
     81   G4cout << "*******************************************************************************" << G4endl;
     82   G4cout << G4endl;
     83
    7384}
    7485

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungContinuous.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeBremsstrahlungContinuous.cc,v 1.10 2008/12/09 15:08:13 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeBremsstrahlungContinuous.cc,v 1.11 2008/12/15 09:23:06 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    9797 G4String dirFile = pathString + "/penelope/" + name;
    9898 std::ifstream file(dirFile);
    99  std::filebuf* lsdp = file.rdbuf();
    100  if (!(lsdp->is_open()))
     99 if (!file.is_open())
    101100     {
    102101      G4String excep = "G4PenelopeBremsstrahlungContinuous - data file " + name + " not found!";
    … …  
    112111     for (size_t j=0;j<NumberofKPoints;j++){
    113112       file >> a1;
    114        ReducedCS[i][j]=a1*cm2; //coversion present in Penelope source
     113       ReducedCS[i][j]=a1*cm2;
    115114     }
    116115     //3) read the total cross section, in cm2
    117116     file >> a1;
    118      TotalCS[i]=a1*cm2; //conversion present in Penelope source
     117     TotalCS[i]=a1*cm2;
    119118     // Check closing item
    120119     file >> a1;
    121120     if (a1 != ((G4double) -1))
    122121       {
    123          G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file " + name;
     122         G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file "
     123           + name;
    124124         G4Exception(excep);
    125125       }
    … …  
    229229 
    230230  //Global x-section factor
    231   G4double Fact=Zmat*Zmat*(energy+electron_mass_c2)*(energy+electron_mass_c2)/(energy*(energy+2.0*electron_mass_c2));
     231  G4double Fact=Zmat*Zmat*(energy+electron_mass_c2)*(energy+electron_mass_c2)/
     232    (energy*(energy+2.0*electron_mass_c2));
    232233  Fact *= PositronCorrection(energy);
    233234

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeCompton.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeCompton.cc,v 1.33 2008/06/03 15:44:25 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeCompton.cc,v 1.34 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    120120             << G4endl;
    121121    }
     122
     123
     124   G4cout << G4endl;
     125   G4cout << "*******************************************************************************" << G4endl;
     126   G4cout << "*******************************************************************************" << G4endl;
     127   G4cout << "   The class G4PenelopeCompton is NOT SUPPORTED ANYMORE. " << G4endl;
     128   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     129   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     130   G4cout << "*******************************************************************************" << G4endl;
     131   G4cout << "*******************************************************************************" << G4endl;
     132   G4cout << G4endl;
     133
    122134}
    123135

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeComptonModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeComptonModel.cc,v 1.2 2008/12/04 14:11:21 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeComptonModel.cc,v 1.4 2009/04/18 18:29:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    4040//                            that are read from the external files.
    4141// 24 Nov 2008   L Pandola    Find a cleaner way to delete vectors.
     42// 16 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     43//                  - apply internal high-energy limit only in constructor
     44//                  - do not apply low-energy limit (default is 0)
     45//                  - do not apply production threshold on level of the model
    4246//
    4347
    … …  
    6973  fIntrinsicLowEnergyLimit = 100.0*eV;
    7074  fIntrinsicHighEnergyLimit = 100.0*GeV;
    71   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     75  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    7276  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    7377  //
    … …  
    118122
    119123void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* particle,
    120                                        const G4DataVector& cuts)
     124                                        const G4DataVector& cuts)
    121125{
    122126  if (verboseLevel > 3)
    … …  
    124128
    125129  InitialiseElementSelectors(particle,cuts);
    126   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    127     {
    128       G4cout << "G4PenelopeComptonModel: low energy limit increased from " <<
    129         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
    130       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    131     }
    132   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    133     {
    134       G4cout << "G4PenelopeComptonModel: high energy limit decreased from " <<
    135         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
    136       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    137     }
    138 
    139   G4cout << "Penelope Compton model is initialized " << G4endl
    140          << "Energy range: "
    141          << LowEnergyLimit() / keV << " keV - "
    142          << HighEnergyLimit() / GeV << " GeV"
    143          << G4endl;
     130
     131  if (verboseLevel > 0) {
     132    G4cout << "Penelope Compton model is initialized " << G4endl
     133           << "Energy range: "
     134           << LowEnergyLimit() / keV << " keV - "
     135           << HighEnergyLimit() / GeV << " GeV"
     136           << G4endl;
     137  }
    144138
    145139  if(isInitialised) return;
    146 
    147   if(pParticleChange)
    148     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    149   else
    150     fParticleChange = new G4ParticleChangeForGamma();
     140  fParticleChange = GetParticleChangeForGamma();
    151141  isInitialised = true;
    152142}
    … …  
    249239  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
    250240
    251   if (photonEnergy0 <= LowEnergyLimit())
    252   {
     241  if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
     242    {
    253243      fParticleChange->ProposeTrackStatus(fStopAndKill);
    254244      fParticleChange->SetProposedKineticEnergy(0.);
    255245      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
    256246      return ;
    257   }
     247    }
    258248
    259249  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
    … …  
    263253    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
    264254  // atom can be selected effitiantly if element selectors are initialised
    265   const G4Element* anElement = SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy0);
     255  const G4Element* anElement =
     256    SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy0);
    266257  G4int Z = (G4int) anElement->GetZ();
    267258  if (verboseLevel > 2)
    … …  
    320311      }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
    321312    }
    322   else //photonEnergy0<5 MeV
     313  else //photonEnergy0 < 5 MeV
    323314    {
    324315      //Incoherent scattering function for theta=PI
    … …  
    472463  }
    473464 
    474 
    475465  //Create scattered electron
    476466  G4double diffEnergy = photonEnergy0*(1-epsilon);
    477467  ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
    478   G4double Q2 = photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
     468  G4double Q2 =
     469    photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
    479470  G4double cosThetaE; //scattering angle for the electron
    480471  if (Q2 > 1.0e-12)
    … …  
    496487  G4int shellId = shell->ShellId();
    497488  G4double ionEnergyInPenelopeDatabase = ionEnergy;
    498   ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase); //protection against energy non-conservation
    499 
    500   G4double eKineticEnergy = diffEnergy - ionEnergy; //subtract the excitation energy. If not emitted by fluorescence,
     489  //protection against energy non-conservation
     490  ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase); 
     491
     492  //subtract the excitation energy. If not emitted by fluorescence
    501493  //the ionization energy is deposited as local energy deposition
     494  G4double eKineticEnergy = diffEnergy - ionEnergy;
    502495  G4double localEnergyDeposit = ionEnergy;
    503496  G4double energyInFluorescence = 0.; //testing purposes only
    … …  
    505498  if (eKineticEnergy < 0)
    506499    {
    507       //It means that there was some problem/mismatch between the two databases. Try to make it work
     500      //It means that there was some problem/mismatch between the two databases.
     501      //Try to make it work
    508502      //In this case available Energy (diffEnergy) < ionEnergy
    509503      //Full residual energy is deposited locally
    … …  
    511505      eKineticEnergy = 0.0;
    512506    }
    513   G4double cutForLowEnergySecondaryPhotons = 250.0*eV;
    514 
    515   //Get the cut for electrons
    516   const G4ProductionCutsTable* theCoupleTable=
    517     G4ProductionCutsTable::GetProductionCutsTable();
    518   size_t indx = couple->GetIndex();
    519   G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
    520   cutE = std::max(cutForLowEnergySecondaryPhotons,cutE);
    521507     
    522508  //the local energy deposit is what remains: part of this may be spent for fluorescence.
    523   if (fUseAtomicDeexcitation)
    524     {
    525       G4int nPhotons=0;
    526 
    527       G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
    528       cutg = std::max(cutForLowEnergySecondaryPhotons,cutg);
    529 
    530       G4DynamicParticle* aPhoton;
    531       G4AtomicDeexcitation deexcitationManager;
     509  if(DeexcitationFlag() && Z > 5) {
     510
     511    const G4ProductionCutsTable* theCoupleTable=
     512      G4ProductionCutsTable::GetProductionCutsTable();
     513
     514    size_t index = couple->GetIndex();
     515    G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
     516    G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
     517
     518    // Generation of fluorescence
     519    // Data in EADL are available only for Z > 5
     520    // Protection to avoid generating photons in the unphysical case of
     521    // shell binding energy > photon energy
     522    if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
     523      {
     524        G4DynamicParticle* aPhoton;
     525        G4AtomicDeexcitation deexcitationManager;
     526        deexcitationManager.SetCutForSecondaryPhotons(cutg);
     527        deexcitationManager.SetCutForAugerElectrons(cute);
     528        deexcitationManager.ActivateAugerElectronProduction(false);
    532529     
    533       if (Z>5 && (localEnergyDeposit > cutg || localEnergyDeposit > cutE))
    534         {
    535           photonVector = deexcitationManager.GenerateParticles(Z,shellId);
    536           for (size_t k=0;k<photonVector->size();k++){
    537             aPhoton = (*photonVector)[k];
    538             if (aPhoton)
     530        photonVector = deexcitationManager.GenerateParticles(Z,shellId);
     531        if(photonVector)
     532          {
     533            size_t nPhotons = photonVector->size();
     534            for (size_t k=0; k<nPhotons; k++)
    539535              {
    540                 G4double itsCut = cutg;
    541                 if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
    542                 G4double itsEnergy = aPhoton->GetKineticEnergy();
    543                 if (itsEnergy > itsCut && itsEnergy <= localEnergyDeposit)
     536                aPhoton = (*photonVector)[k];
     537                if (aPhoton)
    544538                  {
    545                     nPhotons++;
    546                     localEnergyDeposit -= itsEnergy;
    547                     energyInFluorescence += itsEnergy;
    548                   }
    549                 else
    550                   {
    551                     delete aPhoton;
    552                     (*photonVector)[k]=0;
     539                    G4double itsEnergy = aPhoton->GetKineticEnergy();
     540                    if (itsEnergy <= localEnergyDeposit)
     541                      {
     542                        localEnergyDeposit -= itsEnergy;
     543                        if (aPhoton->GetDefinition() == G4Gamma::Gamma())
     544                          energyInFluorescence += itsEnergy;;
     545                        fvect->push_back(aPhoton);                 
     546                      }
     547                    else
     548                      {
     549                        delete aPhoton;
     550                        (*photonVector)[k]=0;
     551                      }
    553552                  }
    554553              }
    555           }
    556         }
    557     }
     554            delete photonVector;
     555          }
     556      }
     557  }
    558558
    559559  //Produce explicitely the electron only if its proposed kinetic energy is
    560560  //above the cut, otherwise add local energy deposition
    561561  G4DynamicParticle* electron = 0;
    562   if (eKineticEnergy > cutE)
     562  //  if (eKineticEnergy > cutE) // VI: may be consistency problem if cut is applied here
     563  if (eKineticEnergy > 0.0)
    563564    {
    564565      G4double xEl = sinThetaE * std::cos(phi+pi);
    … …  
    576577    }
    577578
    578   //This block below is executed only if there is at least one secondary photon produced by
    579   //AtomicDeexcitation
    580   if (photonVector)
    581     {
    582       for (size_t ll=0;ll<photonVector->size();ll++)
    583         {
    584           if ((*photonVector)[ll])
    585             {
    586               G4DynamicParticle* aFluorescencePhoton = (*photonVector)[ll];
    587               fvect->push_back(aFluorescencePhoton);
    588             }
    589         }
    590     }
    591   delete photonVector;
    592579  if (localEnergyDeposit < 0)
    593580    {

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversion.cc

    r819 r1055  
    9797              << G4endl;
    9898     }
     99
     100   G4cout << G4endl;
     101   G4cout << "*******************************************************************************" << G4endl;
     102   G4cout << "*******************************************************************************" << G4endl;
     103   G4cout << "   The class G4PenelopeGammaConversion is NOT SUPPORTED ANYMORE. " << G4endl;
     104   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     105   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     106   G4cout << "*******************************************************************************" << G4endl;
     107   G4cout << "*******************************************************************************" << G4endl;
     108   G4cout << G4endl;
    99109}
    100110 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversionModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeGammaConversionModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeGammaConversionModel.cc,v 1.4 2009/05/19 14:57:01 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    3232// --------
    3333// 06 Oct 2008   L Pandola    Migration from process to model
     34// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     35//                  - apply internal high-energy limit only in constructor
     36//                  - do not apply low-energy limit (default is 0)
     37//                  - do not apply production threshold on level of the model
     38// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
     39//                            Initialise(), since they might be checked later on
    3440//
    3541
    … …  
    5662  fSmallEnergy = 1.1*MeV;
    5763
    58   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     64  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    5965  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    6066  //
    … …  
    8995      crossSectionHandler->Clear();
    9096      delete crossSectionHandler;
     97      crossSectionHandler = 0;
    9198    }
    9299 
    93   //Check energy limits
    94   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    95     {
    96       G4cout << "G4PenelopeGammaConversionModel: low energy limit increased from " <<
    97         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
    98       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    99     }
    100   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    101     {
    102       G4cout << "G4PenelopeGammaConversionModel: high energy limit decreased from " <<
    103         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
    104       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    105     }
    106 
    107100  //Re-initialize cross section handler
    108101  crossSectionHandler = new G4CrossSectionHandler();
    109   crossSectionHandler->Initialise(0,LowEnergyLimit(),HighEnergyLimit(),400);
     102  crossSectionHandler->Initialise(0,fIntrinsicLowEnergyLimit,HighEnergyLimit(),400);
    110103  crossSectionHandler->Clear();
    111104  G4String crossSectionFile = "penelope/pp-cs-pen-";
    … …  
    117110    G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
    118111
    119   G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
    120          << "Energy range: "
    121          << LowEnergyLimit() / MeV << " MeV - "
    122          << HighEnergyLimit() / GeV << " GeV"
    123          << G4endl;
     112  if (verboseLevel > 0) {
     113    G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
     114           << "Energy range: "
     115           << LowEnergyLimit() / MeV << " MeV - "
     116           << HighEnergyLimit() / GeV << " GeV"
     117           << G4endl;
     118  }
    124119
    125120  if(isInitialised) return;
    126 
    127   if(pParticleChange)
    128     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    129   else
    130     fParticleChange = new G4ParticleChangeForGamma();
     121  fParticleChange = GetParticleChangeForGamma();
    131122  isInitialised = true;
    132123}
    … …  
    151142
    152143  G4int iZ = (G4int) Z;
    153   if (!crossSectionHandler)
    154     {
    155       G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl;
    156       G4cout << "The cross section handler is not correctly initialized" << G4endl;
    157       G4Exception();
    158     }
     144  //  if (!crossSectionHandler) //VI: should not be checked in run time
     145  //  {
     146  //    G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl;
     147  //    G4cout << "The cross section handler is not correctly initialized" << G4endl;
     148  //    G4Exception();
     149  //  }
    159150  G4double cs = crossSectionHandler->FindValue(iZ,energy);
    160151
    … …  
    167158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    168159
    169 void G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
    170                                               const G4MaterialCutsCouple* couple,
    171                                               const G4DynamicParticle* aDynamicGamma,
    172                                               G4double,
    173                                               G4double)
     160void
     161G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
     162                                                  const G4MaterialCutsCouple* couple,
     163                                                  const G4DynamicParticle* aDynamicGamma,
     164                                                  G4double,
     165                                                  G4double)
    174166{
    175167  //
    … …  
    191183  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
    192184
    193   if (photonEnergy <= LowEnergyLimit())
    194   {
    195       fParticleChange->ProposeTrackStatus(fStopAndKill);
    196       fParticleChange->SetProposedKineticEnergy(0.);
     185  // Always kill primary
     186  fParticleChange->ProposeTrackStatus(fStopAndKill);
     187  fParticleChange->SetProposedKineticEnergy(0.);
     188
     189  if (photonEnergy <= fIntrinsicLowEnergyLimit)
     190    {
    197191      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
    198192      return ;
    199   }
     193    }
    200194
    201195  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
    … …  
    310304
    311305  //Generate explicitely the electron in the pair, only if it is > threshold
    312   const G4ProductionCutsTable* theCoupleTable=
    313     G4ProductionCutsTable::GetProductionCutsTable();
    314   size_t indx = couple->GetIndex();
    315   G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
    316   //G4double cutP = (*(theCoupleTable->GetEnergyCutsVector(2)))[indx];
    317 
    318   if (electronKineEnergy > cutE)
     306  //VI: applying cut here provides inconsistency
     307
     308  if (electronKineEnergy > 0.0)
    319309    {
    320310      G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
    … …  
    333323  //Generate the positron. Real particle in any case, because it will annihilate. If below
    334324  //threshold, produce it at rest
    335   if (positronKineEnergy < cutE)
     325  // VI: here there was a bug - positron and electron cuts are different
     326  if (positronKineEnergy < 0.0)
    336327    {
    337328      localEnergyDeposit += positronKineEnergy;
    … …  
    344335  fvect->push_back(positron);
    345336
    346   //Update the status of the primary gamma (kill it)
    347   fParticleChange->SetProposedKineticEnergy(0.);
     337  //Add rest of energy to the local energy deposit
    348338  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
    349   fParticleChange->ProposeTrackStatus(fStopAndKill);
    350 
    351339 
    352  if (verboseLevel > 1)
     340  if (verboseLevel > 1)
    353341    {
    354342      G4cout << "-----------------------------------------------------------" << G4endl;
    … …  
    357345      G4cout << "-----------------------------------------------------------" << G4endl;
    358346      if (electronKineEnergy)
    359         G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" << G4endl;
     347        G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
     348               << G4endl;
    360349      if (positronKineEnergy)
    361350        G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
    … …  
    373362                                      localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
    374363      if (energyDiff > 0.05*keV)
    375         G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " <<
    376           (electronKineEnergy+positronKineEnergy+
    377            localEnergyDeposit+2.0*electron_mass_c2)/keV << " keV (final) vs. " <<
    378           photonEnergy/keV << " keV (initial)" << G4endl;
    379     }
    380  
     364        G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: "
     365               << (electronKineEnergy+positronKineEnergy+
     366                   localEnergyDeposit+2.0*electron_mass_c2)/keV
     367               << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
     368    }
    381369}
    382370

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeIonisation.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeIonisation.cc,v 1.19 2006/06/29 19:40:49 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeIonisation.cc,v 1.20 2009/05/02 09:59:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// --------------------------------------------------------------
    … …  
    9393  ReadData(); //Read data from file
    9494 
     95   G4cout << G4endl;
     96   G4cout << "*******************************************************************************" << G4endl;
     97   G4cout << "*******************************************************************************" << G4endl;
     98   G4cout << "   The class G4PenelopeIonisation is NOT SUPPORTED ANYMORE. " << G4endl;
     99   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     100   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     101   G4cout << "*******************************************************************************" << G4endl;
     102   G4cout << "*******************************************************************************" << G4endl;
     103   G4cout << G4endl;
     104
    95105}
    96106

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeIonisationModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeIonisationModel.cc,v 1.2 2008/12/05 09:15:43 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeIonisationModel.cc,v 1.5 2009/05/19 14:57:01 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    3232// --------
    3333// 26 Nov 2008   L Pandola    Migration from process to model
     34// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     35//                  - apply internal high-energy limit only in constructor
     36//                  - do not apply low-energy limit (default is 0)
     37//                  - added MinEnergyCut method
     38//                  - do not change track status
     39// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
     40//                            Initialise(), since they might be checked later on
    3441//
    3542
    … …  
    4855#include "G4CrossSectionHandler.hh"
    4956#include "G4AtomicDeexcitation.hh"
    50 #include "G4ProcessManager.hh"
    5157#include "G4VEMDataSet.hh"
    5258
    … …  
    6874  fIntrinsicLowEnergyLimit = 100.0*eV;
    6975  fIntrinsicHighEnergyLimit = 100.0*GeV;
    70   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     76  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    7177  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    7278  //
    … …  
    136142    G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
    137143
    138   InitialiseElementSelectors(particle,cuts);
    139 
    140144  //Delete and re-initialize the cross section handler
    141145  if (crossSectionHandler)
    … …  
    143147      crossSectionHandler->Clear();
    144148      delete crossSectionHandler;
     149      crossSectionHandler = 0;
    145150    }
    146151
    … …  
    148153    {
    149154      for (size_t i=0; i<theXSTable->size(); i++)
    150         delete (*theXSTable)[i];
     155        {
     156          delete (*theXSTable)[i];
     157          (*theXSTable)[i] = 0;
     158        }
    151159      delete theXSTable;
     160      theXSTable = 0;
    152161     }
    153 
    154 
    155   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    156     {
    157       G4cout << "G4PenelopeIonisationModel: low energy limit increased from " <<
    158         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
    159       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    160     }
    161   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    162     {
    163       G4cout << "G4PenelopeIonisationModel: high energy limit decreased from " <<
    164         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
    165       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    166     }
    167162
    168163  crossSectionHandler = new G4CrossSectionHandler();
    … …  
    177172  //This is used to retrieve cross section values later on
    178173  crossSectionHandler->BuildMeanFreePathForMaterials();
    179  
    180  if (verboseLevel > 2)
     174
     175  InitialiseElementSelectors(particle,cuts);
     176 
     177  if (verboseLevel > 2)
    181178    G4cout << "Loaded cross section files for PenelopeIonisationModel" << G4endl;
    182179
    183 
    184   G4cout << "Penelope Ionisation model is initialized " << G4endl
    185          << "Energy range: "
    186          << LowEnergyLimit() / keV << " keV - "
    187          << HighEnergyLimit() / GeV << " GeV"
    188          << G4endl;
     180  if (verboseLevel > 2) {
     181    G4cout << "Penelope Ionisation model is initialized " << G4endl
     182           << "Energy range: "
     183           << LowEnergyLimit() / keV << " keV - "
     184           << HighEnergyLimit() / GeV << " GeV"
     185           << G4endl;
     186  }
    189187
    190188  if(isInitialised) return;
    191 
    192   if(pParticleChange)
    193     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
    194   else
    195     fParticleChange = new G4ParticleChangeForLoss();
     189  fParticleChange = GetParticleChangeForLoss();
    196190  isInitialised = true;
    197191}
    … …  
    227221  SetupForMaterial(theParticle, material, energy);
    228222
     223  // VI - should be check at initialisation not in run time
     224  /*
    229225  if (!crossSectionHandler)
    230226    {
    … …  
    233229      G4Exception();
    234230    }
    235  
     231  */ 
    236232  if (!theXSTable)
    237233    {
    … …  
    244240      theXSTable = BuildCrossSectionTable(theParticle);
    245241    }
    246 
     242 
    247243  G4double totalCross = 0.0;
    248244  G4double cross = 0.0;
    … …  
    356352                                              const G4MaterialCutsCouple* couple,
    357353                                              const G4DynamicParticle* aDynamicParticle,
    358                                               G4double,
    359                                               G4double)
     354                                              G4double cutE, G4double)
    360355{
    361356  // Penelope model to sample the final state following an hard inelastic interaction.
    … …  
    397392  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
    398393
    399   if (kineticEnergy0 <= LowEnergyLimit())
     394  if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
    400395  {
    401396    fParticleChange->SetProposedKineticEnergy(0.);
    402397    fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
    403     //Check if there are AtRest processes
    404     if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
    405       //In this case there is at least one AtRest process
    406       fParticleChange->ProposeTrackStatus(fStopButAlive);
    407     else
    408       fParticleChange->ProposeTrackStatus(fStopAndKill);
    409398    return ;
    410399  }
    … …  
    430419  G4int iZ = SampleRandomAtom(couple,kineticEnergy0);
    431420 
    432   G4double cutForLowEnergySecondaryParticles = 250.0*eV;
    433421  const G4ProductionCutsTable* theCoupleTable=
    434422    G4ProductionCutsTable::GetProductionCutsTable();
    435423  size_t indx = couple->GetIndex();
    436   G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
    437424  G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
    438 
    439   //Production cut for delta-rays (electrons)
    440   cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
    441   //Production cut for gamma (fluorescence)
    442   cutG = std::max(cutForLowEnergySecondaryParticles,cutG);
    443425 
    444426  if (verboseLevel > 2)
    … …  
    474456  else
    475457    {
    476       fParticleChange->ProposeMomentumDirection(electronDirection1);
    477       fParticleChange->SetProposedKineticEnergy(0.*eV);
    478       if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
    479         //In this case there is at least one AtRest process
    480         fParticleChange->ProposeTrackStatus(fStopButAlive);
    481       else
    482         fParticleChange->ProposeTrackStatus(fStopAndKill);
     458      fParticleChange->SetProposedKineticEnergy(0.);
    483459    }
    484460 
    … …  
    517493    //Penelope subtracted the fluorescence, but one has to match the databases
    518494    eKineticEnergy = energySecondary+ioniEnergy-bindingEnergy;
    519    
    520   //VERIFICA QUI LA STORIA DEL LOCAL ENERGY DEPOSIT!
    521 
     495 
    522496  G4double localEnergyDeposit = ionEnergy;
    523497  G4double energyInFluorescence = 0.0*eV;
    524498
    525    std::vector<G4DynamicParticle*> *photonVector = 0;
    526   if (fUseAtomicDeexcitation)
    527     {
    528       if (iZ>5 && (ionEnergy > cutG || ionEnergy > cutE))
     499  if(DeexcitationFlag() && iZ > 5)
     500    {
     501      if (ionEnergy > cutG || ionEnergy > cutE)
    529502        {
    530           photonVector = deexcitationManager.GenerateParticles(iZ,shellId);
    531           //Check for single photons if they are above threshold
    532           for (size_t k=0;k<photonVector->size();k++)
     503          deexcitationManager.SetCutForSecondaryPhotons(cutG);
     504          deexcitationManager.SetCutForAugerElectrons(cutE);
     505          std::vector<G4DynamicParticle*> *photonVector =
     506            deexcitationManager.GenerateParticles(iZ,shellId);
     507          //Check for secondaries
     508          if(photonVector)
    533509            {
    534               G4DynamicParticle* aPhoton = (*photonVector)[k];
    535               if (aPhoton)
     510              for (size_t k=0;k<photonVector->size();k++)
    536511                {
    537                   G4double itsCut = cutG;
    538                   if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
    539                   G4double itsEnergy = aPhoton->GetKineticEnergy();
    540                   if (itsEnergy > itsCut && itsEnergy <= ionEnergy)
     512                  G4DynamicParticle* aPhoton = (*photonVector)[k];
     513                  if (aPhoton)
    541514                    {
    542                       localEnergyDeposit -= itsEnergy;
    543                       energyInFluorescence += itsEnergy;
    544                     }
    545                   else
    546                     {
    547                       delete aPhoton;
    548                       (*photonVector)[k] = 0;
     515                      G4double itsEnergy = aPhoton->GetKineticEnergy();
     516                      if (itsEnergy <= localEnergyDeposit)
     517                        {
     518                          if(aPhoton->GetDefinition() == G4Gamma::Gamma())
     519                            energyInFluorescence += itsEnergy;
     520                          localEnergyDeposit -= itsEnergy;
     521                          fvect->push_back(aPhoton);
     522                        }
     523                      else
     524                        {
     525                          delete aPhoton;
     526                          (*photonVector)[k] = 0;
     527                        }
    549528                    }
    550529                }
     530              delete photonVector;
    551531            }
    552532        }
    … …  
    567547  fvect->push_back(deltaElectron);
    568548 
    569   //Generate fluorescence, if it is the case
    570   //This block is executed only if there is at least one secondary photon produced by
    571   //G4AtomicDeexcitation
    572   if (photonVector)
    573     {
    574       for (size_t ll=0;ll<photonVector->size();ll++)
    575         if ((*photonVector)[ll])
    576           {
    577             G4DynamicParticle* aFluorescencePhoton = (*photonVector)[ll];
    578             fvect->push_back(aFluorescencePhoton);
    579           }
    580     }
    581   delete photonVector;
    582 
    583549  if (localEnergyDeposit < 0)
    584550    {
    … …  
    753719//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    754720
    755 void G4PenelopeIonisationModel::CalculateDiscreteForElectrons(G4double kinEnergy,G4double cutoffEnergy,
    756                                                               G4int Z,G4double electronVolumeDensity)
     721void
     722G4PenelopeIonisationModel::CalculateDiscreteForElectrons(G4double kinEnergy,
     723                                                         G4double cutoffEnergy,
     724                                                         G4int Z,
     725                                                         G4double electronVolumeDensity)
    757726{
    758727  if (verboseLevel > 2)
    … …  
    1030999//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    10311000
    1032  void G4PenelopeIonisationModel::CalculateDiscreteForPositrons(G4double kinEnergy,G4double cutoffEnergy,
    1033                                                                G4int Z,G4double electronVolumeDensity)
     1001 void
     1002G4PenelopeIonisationModel::CalculateDiscreteForPositrons(G4double kinEnergy,
     1003                                                         G4double cutoffEnergy,
     1004                                                         G4int Z,
     1005                                                         G4double electronVolumeDensity)
    10341006{
    10351007  kineticEnergy1=kinEnergy;
    … …  
    12871259//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    12881260
    1289 G4double G4PenelopeIonisationModel::CalculateCrossSectionsRatio(G4double kinEnergy,
    1290                                                                 G4double cutoffEnergy,
    1291                                                                 G4int Z, G4double electronVolumeDensity,
    1292                                                                 const G4ParticleDefinition* theParticle)
     1261G4double
     1262G4PenelopeIonisationModel::CalculateCrossSectionsRatio(G4double kinEnergy,
     1263                                                       G4double cutoffEnergy,
     1264                                                       G4int Z, G4double electronVolumeDensity,
     1265                                                       const G4ParticleDefinition* theParticle)
    12931266{
    12941267  //Constants
    … …  
    13201293//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    13211294
    1322 std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForElectrons(G4double kineticEnergy,
    1323                                                                                        G4double resEnergy,
    1324                                                                                        G4double densityCorrection,
    1325                                                                                        G4double cutoffEnergy)
     1295std::pair<G4double,G4double>
     1296G4PenelopeIonisationModel::CrossSectionsRatioForElectrons(G4double kineticEnergy,
     1297                                                          G4double resEnergy,
     1298                                                          G4double densityCorrection,
     1299                                                          G4double cutoffEnergy)
    13261300{
    13271301  std::pair<G4double,G4double> theResult(0.,0.);
    … …  
    14101384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    14111385
    1412 std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForPositrons(G4double kineticEnergy,
    1413                                                                                        G4double resEnergy,
    1414                                                                                        G4double densityCorrection,
    1415                                                                                        G4double cutoffEnergy)
     1386std::pair<G4double,G4double>
     1387G4PenelopeIonisationModel::CrossSectionsRatioForPositrons(G4double kineticEnergy,
     1388                                                          G4double resEnergy,
     1389                                                          G4double densityCorrection,
     1390                                                          G4double cutoffEnergy)
    14161391{
    14171392
    … …  
    16401615#include "G4EMDataSet.hh"
    16411616
    1642 std::vector<G4VEMDataSet*>* G4PenelopeIonisationModel::BuildCrossSectionTable(const
    1643                                                                               G4ParticleDefinition* theParticle)
     1617std::vector<G4VEMDataSet*>*
     1618G4PenelopeIonisationModel::BuildCrossSectionTable(const G4ParticleDefinition* theParticle)
    16441619{
    16451620  std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectric.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4PenelopePhotoElectric.cc,v 1.12 2006/06/29 19:40:51 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4PenelopePhotoElectric.cc,v 1.14 2009/05/02 09:59:17 sincerti Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Author: L. Pandola
    … …  
    3737// 31 May 2005  L. Pandola  Added Sauter formula for the sampling of
    3838//                          the electron direction
     39// 08 Jan 2009  L. Pandola  Check shell index to avoid mismatch between
     40//                          the Penelope cross section database and the
     41//                          G4AtomicTransitionManager database. It suppresses
     42//                          a warning from G4AtomicTransitionManager only.
     43//                          Results are unchanged.
    3944// --------------------------------------------------------------
    4045
    … …  
    8994             << G4endl;
    9095    }
     96
     97   G4cout << G4endl;
     98   G4cout << "*******************************************************************************" << G4endl;
     99   G4cout << "*******************************************************************************" << G4endl;
     100   G4cout << "   The class G4PenelopePhotoElectric is NOT SUPPORTED ANYMORE. " << G4endl;
     101   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     102   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     103   G4cout << "*******************************************************************************" << G4endl;
     104   G4cout << "*******************************************************************************" << G4endl;
     105   G4cout << G4endl;
     106
    91107}
    92108
    … …  
    146162  // Retrieve the corresponding identifier and binding energy of the selected shell
    147163  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
     164
     165  //The number of shell cross section possibly reported in the Penelope database
     166  //might be different from the number of shells in the G4AtomicTransitionManager
     167  //(namely, Penelope may contain more shell, especially for very light elements).
     168  //In order to avoid a warning message from the G4AtomicTransitionManager, I
     169  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
     170  //has a shellID>maxID, it sets the shellID to the last valid shell.
     171  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
     172  if (shellIndex >= numberOfShells)
     173    shellIndex = numberOfShells-1;
     174
    148175  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
    149176  G4double bindingEnergy = shell->BindingEnergy();

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectricModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopePhotoElectricModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopePhotoElectricModel.cc,v 1.6 2009/05/19 14:57:01 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    3131// History:
    3232// --------
    33 // 08 Oct 2008   L Pandola    Migration from process to model
     33// 08 Oct 2008   L Pandola  Migration from process to model
     34// 08 Jan 2009   L Pandola  Check shell index to avoid mismatch between
     35//                          the Penelope cross section database and the
     36//                          G4AtomicTransitionManager database. It suppresses
     37//                          a warning from G4AtomicTransitionManager only.
     38//                          Results are unchanged.
     39// 25 Mar 2009   L Pandola  Small fix to avoid wrong energy-violation warnings
     40// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     41//                  - apply internal high-energy limit only in constructor
     42//                  - do not apply low-energy limit (default is 0)
     43//                  - do not apply production threshold on secondaries
     44// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
     45//                            Initialise(), since they might be checked later on
    3446//
    3547
    … …  
    5870  fIntrinsicLowEnergyLimit = 100.0*eV;
    5971  fIntrinsicHighEnergyLimit = 100.0*GeV;
    60   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     72  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    6173  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    6274  //
    … …  
    90102      crossSectionHandler->Clear();
    91103      delete crossSectionHandler;
     104      crossSectionHandler = 0;
    92105    }
    93106  if (shellCrossSectionHandler)
    … …  
    95108      shellCrossSectionHandler->Clear();
    96109      delete shellCrossSectionHandler;
    97     }
    98 
    99   //Check energy limits
    100   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    101     {
    102       G4cout << "G4PenelopePhotoElectricModel: low energy limit increased from " <<
    103         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" <<
    104         G4endl;
    105       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    106     }
    107   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    108     {
    109       G4cout << "G4PenelopePhotoElectricModel: high energy limit decreased from " <<
    110         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV"
    111              << G4endl;
    112       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    113     }
    114 
     110      shellCrossSectionHandler =0;
     111    }
    115112
    116113  //Re-initialize cross section handlers
    … …  
    129126    G4cout << "Loaded cross section files for PenelopePhotoElectric" << G4endl;
    130127
    131   G4cout << "Penelope Photo-Electric model is initialized " << G4endl
    132          << "Energy range: "
    133          << LowEnergyLimit() / MeV << " MeV - "
    134          << HighEnergyLimit() / GeV << " GeV"
    135          << G4endl;
     128  if (verboseLevel > 0) {
     129    G4cout << "Penelope Photo-Electric model is initialized " << G4endl
     130           << "Energy range: "
     131           << LowEnergyLimit() / MeV << " MeV - "
     132           << HighEnergyLimit() / GeV << " GeV"
     133           << G4endl;
     134  }
    136135
    137136  if(isInitialised) return;
    138 
    139   if(pParticleChange)
    140     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    141   else
    142     fParticleChange = new G4ParticleChangeForGamma();
     137  fParticleChange = GetParticleChangeForGamma();
    143138  isInitialised = true;
    144139}
    … …  
    162157
    163158  G4int iZ = (G4int) Z;
    164   if (!crossSectionHandler)
    165     {
    166       G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
    167       G4cout << "The cross section handler is not correctly initialized" << G4endl;
    168       G4Exception();
    169     }
     159  //  if (!crossSectionHandler) // VI: should not be
     160  //  {
     161  //    G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
     162  //    G4cout << "The cross section handler is not correctly initialized" << G4endl;
     163  //    G4Exception();
     164  //  }
    170165  G4double cs = crossSectionHandler->FindValue(iZ,energy);
    171166 
    … …  
    203198  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
    204199
    205   if (photonEnergy <= LowEnergyLimit())
    206   {
    207       fParticleChange->ProposeTrackStatus(fStopAndKill);
    208       fParticleChange->SetProposedKineticEnergy(0.);
     200  // always kill primary
     201  fParticleChange->ProposeTrackStatus(fStopAndKill);
     202  fParticleChange->SetProposedKineticEnergy(0.);
     203
     204  if (photonEnergy <= fIntrinsicLowEnergyLimit)
     205    {
    209206      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
    210207      return ;
    211   }
     208    }
    212209
    213210  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
    … …  
    217214    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
    218215  //use crossSectionHandler instead of G4EmElementSelector because in this case
    219   //the dimension of the table is equal to the dimension of the database (less interpolation errors)
     216  //the dimension of the table is equal to the dimension of the database
     217  //(less interpolation errors)
    220218  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
    221219  if (verboseLevel > 2)
    … …  
    227225  // Retrieve the corresponding identifier and binding energy of the selected shell
    228226  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
     227
     228  //The number of shell cross section possibly reported in the Penelope database
     229  //might be different from the number of shells in the G4AtomicTransitionManager
     230  //(namely, Penelope may contain more shell, especially for very light elements).
     231  //In order to avoid a warning message from the G4AtomicTransitionManager, I
     232  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
     233  //has a shellID>maxID, it sets the shellID to the last valid shell.
     234  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
     235  if (shellIndex >= numberOfShells)
     236    shellIndex = numberOfShells-1;
     237
    229238  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
    230239  G4double bindingEnergy = shell->BindingEnergy();
    … …  
    236245  G4double eKineticEnergy = photonEnergy - bindingEnergy;
    237246 
    238   G4double cutForLowEnergySecondaryParticles = 250.0*eV;
    239   const G4ProductionCutsTable* theCoupleTable=
    240     G4ProductionCutsTable::GetProductionCutsTable();
    241   size_t indx = couple->GetIndex();
    242   G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
    243   cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
    244 
    245247  // There may be cases where the binding energy of the selected shell is > photon energy
    246248  // In such cases do not generate secondaries
    … …  
    248250    {
    249251      //Now check if the electron is above cuts: if so, it is created explicitely
    250       if (eKineticEnergy > cutE)
    251         {
    252           // The electron is created
    253           // Direction sampled from the Sauter distribution
    254           G4double cosTheta = SampleElectronDirection(eKineticEnergy);
    255           G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
    256           G4double phi = twopi * G4UniformRand() ;
    257           G4double dirx = sinTheta * std::cos(phi);
    258           G4double diry = sinTheta * std::sin(phi);
    259           G4double dirz = cosTheta ;
    260           G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
    261           electronDirection.rotateUz(photonDirection);
    262           G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
    263                                                                electronDirection,
    264                                                                eKineticEnergy);
    265           fvect->push_back(electron);
    266         }
    267       else
    268         localEnergyDeposit += eKineticEnergy;   
    269     }
     252      //VI: checking cut here provides inconsistency in testing
     253      //      if (eKineticEnergy > cutE)
     254      // The electron is created
     255      // Direction sampled from the Sauter distribution
     256      G4double cosTheta = SampleElectronDirection(eKineticEnergy);
     257      G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
     258      G4double phi = twopi * G4UniformRand() ;
     259      G4double dirx = sinTheta * std::cos(phi);
     260      G4double diry = sinTheta * std::sin(phi);
     261      G4double dirz = cosTheta ;
     262      G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
     263      electronDirection.rotateUz(photonDirection);
     264      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
     265                                                           electronDirection,
     266                                                           eKineticEnergy);
     267      fvect->push_back(electron);
     268    }
     269  //  else
     270  //  {
     271  //    localEnergyDeposit += eKineticEnergy;   
     272  //    eKineticEnergy = 0;
     273  //  }
     274  //  }
    270275  else
     276    {
    271277      bindingEnergy = photonEnergy;
    272    
    273  
     278    }
    274279  G4double energyInFluorescence = 0; //testing purposes
    275280
    276281  //Now, take care of fluorescence, if required
    277   if (fUseAtomicDeexcitation)
    278     {
     282  if(DeexcitationFlag() && Z > 5)
     283    {
     284      const G4ProductionCutsTable* theCoupleTable=
     285        G4ProductionCutsTable::GetProductionCutsTable();
     286      size_t indx = couple->GetIndex();
    279287      G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
    280       cutG = std::min(cutForLowEnergySecondaryParticles,cutG);
     288      G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
    281289     
    282       std::vector<G4DynamicParticle*>* photonVector = 0;
    283 
    284290      // Protection to avoid generating photons in the unphysical case of
    285291      // shell binding energy > photon energy
    286       if (Z > 5  && (bindingEnergy > cutG || bindingEnergy > cutE))
     292      if (bindingEnergy > cutG || bindingEnergy > cutE)
    287293        {
    288           photonVector = deexcitationManager.GenerateParticles(Z,shellId);
    289           //Check for single photons (if they are above threshold)
    290           for (size_t k=0; k< photonVector->size(); k++)
     294          deexcitationManager.SetCutForSecondaryPhotons(cutG);
     295          deexcitationManager.SetCutForAugerElectrons(cutE);
     296          std::vector<G4DynamicParticle*>* photonVector =
     297            deexcitationManager.GenerateParticles(Z,shellId);
     298          //Check for secondaries
     299          if(photonVector)
    291300            {
    292               G4DynamicParticle* aPhoton = (*photonVector)[k];
    293               if (aPhoton)
     301              for (size_t k=0; k< photonVector->size(); k++)
    294302                {
    295                   G4double itsCut = cutG;
    296                   if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
    297                   G4double itsEnergy = aPhoton->GetKineticEnergy();
    298                  
    299                   if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
     303                  G4DynamicParticle* aPhoton = (*photonVector)[k];
     304                  if (aPhoton)
    300305                    {
    301                       // Local energy deposit is given as the sum of the
    302                       // energies of incident photons minus the energies
    303                       // of the outcoming fluorescence photons
    304                       bindingEnergy -= itsEnergy;
    305                      
    306                     }
    307                   else
    308                     {
    309                       (*photonVector)[k] = 0;
     306                      G4double itsEnergy = aPhoton->GetKineticEnergy();
     307                      if (itsEnergy <= bindingEnergy)
     308                        {
     309                          if(aPhoton->GetDefinition() == G4Gamma::Gamma())
     310                            energyInFluorescence += itsEnergy;
     311                          bindingEnergy -= itsEnergy;
     312                          fvect->push_back(aPhoton);
     313                        }
     314                      else
     315                        {
     316                          delete aPhoton;
     317                          (*photonVector)[k] = 0;
     318                        }
    310319                    }
    311320                }
     321              delete photonVector;
    312322            }
    313         }
    314       //Register valid secondaries
    315       if (photonVector)
    316         {
    317           for ( size_t ll = 0; ll <photonVector->size(); ll++)
    318             {
    319               G4DynamicParticle* aPhoton = (*photonVector)[ll];
    320               if (aPhoton)
    321                 {
    322                   energyInFluorescence += aPhoton->GetKineticEnergy();
    323                   fvect->push_back(aPhoton);
    324                 }
    325             }
    326           delete photonVector;
    327323        }
    328324    }
    … …  
    330326  localEnergyDeposit += bindingEnergy;
    331327     
    332 
    333328  if (localEnergyDeposit < 0)
    334329    {
    … …  
    339334    }
    340335
    341  //Update the status of the primary gamma (kill it)
    342   fParticleChange->SetProposedKineticEnergy(0.);
    343336  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
    344   fParticleChange->ProposeTrackStatus(fStopAndKill);
    345337
    346338  if (verboseLevel > 1)
    … …  
    360352  if (verboseLevel > 0)
    361353    {
    362       G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
     354      G4double energyDiff =
     355        std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
    363356      if (energyDiff > 0.05*keV)
    364357        G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
    365           (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << " keV (final) vs. " <<
     358          (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV
     359               << " keV (final) vs. " <<
    366360          photonEnergy/keV << " keV (initial)" << G4endl;
    367361    }

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeRayleigh.cc

    r1007 r1055  
    2626// --------------------------------------------------------------------
    2727//
    28 // $Id: G4PenelopeRayleigh.cc,v 1.15 2007/09/03 09:43:14 pandola Exp $
    29 // GEANT4 tag $Name: geant4-09-02 $
     28// $Id: G4PenelopeRayleigh.cc,v 1.17 2009/05/02 09:59:17 sincerti Exp $
     29// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    3030//
    3131// Author: L. Pandola (luciano.pandola@cern.ch)
    … …  
    4141// 03 Sep 2007 L.Pandola      Bug fix for the filling of physics table for
    4242//                            compounds defined by the mass fraction (bug #965)
     43// 02 Apr 2009 L.Pandola      Bux fixed in the calculation of mfp for compound
     44//                            materials defined as fraction of mass
     45//                            (reported by Zhang Qiwei)
    4346// --------------------------------------------------------------------
    4447
    … …  
    9396              << G4endl;
    9497     }
     98
     99   G4cout << G4endl;
     100   G4cout << "*******************************************************************************" << G4endl;
     101   G4cout << "*******************************************************************************" << G4endl;
     102   G4cout << "   The class G4PenelopeRayleigh is NOT SUPPORTED ANYMORE. " << G4endl;
     103   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
     104   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
     105   G4cout << "*******************************************************************************" << G4endl;
     106   G4cout << "*******************************************************************************" << G4endl;
     107   G4cout << G4endl;
    95108}
    96109
    … …  
    162175          const G4double* vector_of_atoms = material->GetVecNbOfAtomsPerVolume();
    163176          const G4int* stechiometric = material->GetAtomsVector();
    164           G4double density = vector_of_atoms[iright]; //non-bound molecules (default)
     177          //cs is the cross section _per atom_ in the case of compounds, while it is
     178          //_per molecule_ in the case of molecules
     179          G4double density = material->GetTotNbOfAtomsPerVolume();  //non-bound molecules (default)
    165180          if (stechiometric)
    166181            {

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeRayleighModel.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PenelopeRayleighModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PenelopeRayleighModel.cc,v 1.4 2009/05/19 14:57:01 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    3232// --------
    3333// 14 Oct 2008   L Pandola    Migration from process to model
     34// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
     35//                  - apply internal high-energy limit only in constructor
     36//                  - do not apply low-energy limit (default is 0)
     37// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
     38//                            PrepareConstants(), since they might be checked later on
    3439//
    3540
    … …  
    5661  fIntrinsicLowEnergyLimit = 100.0*eV;
    5762  fIntrinsicHighEnergyLimit = 100.0*GeV;
    58   SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
     63  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    5964  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    6065  //
    … …  
    9095    G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
    9196
    92   if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
    93     {
    94       G4cout << "G4PenelopeRayleighModel: low energy limit increased from " <<
    95         LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
    96       SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
    97     }
    98   if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
    99     {
    100       G4cout << "G4PenelopeRayleighModel: high energy limit decreased from " <<
    101         HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
    102       SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
    103     }
    104 
    105   G4cout << "Penelope Rayleigh model is initialized " << G4endl
    106          << "Energy range: "
    107          << LowEnergyLimit() / keV << " keV - "
    108          << HighEnergyLimit() / GeV << " GeV"
    109          << G4endl;
     97
     98  if (verboseLevel > 0) {
     99    G4cout << "Penelope Rayleigh model is initialized " << G4endl
     100           << "Energy range: "
     101           << LowEnergyLimit() / keV << " keV - "
     102           << HighEnergyLimit() / GeV << " GeV"
     103           << G4endl;
     104  }
    110105
    111106  if(isInitialised) return;
    112 
    113   if(pParticleChange)
    114     fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    115   else
    116     fParticleChange = new G4ParticleChangeForGamma();
     107  fParticleChange = GetParticleChangeForGamma();
    117108  isInitialised = true;
    118109}
    … …  
    120111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    121112
    122 G4double G4PenelopeRayleighModel::CrossSectionPerVolume(const G4Material* material,
    123                                            const G4ParticleDefinition* p,
    124                                            G4double ekin,
    125                                            G4double,
    126                                            G4double)
     113G4double
     114G4PenelopeRayleighModel::CrossSectionPerVolume(const G4Material* material,
     115                                               const G4ParticleDefinition* p,
     116                                               G4double ekin,
     117                                               G4double,
     118                                               G4double)
    127119{
    128120  // Penelope model to calculate the Rayleigh scattering inverse mean
    … …  
    233225
    234226void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
    235                                               const G4MaterialCutsCouple* couple,
    236                                               const G4DynamicParticle* aDynamicGamma,
    237                                               G4double,
    238                                               G4double)
     227                                                const G4MaterialCutsCouple* couple,
     228                                                const G4DynamicParticle* aDynamicGamma,
     229                                                G4double,
     230                                                G4double)
    239231{
    240232  // Penelope model to sample the Rayleigh scattering final state.
    … …  
    249241  if (verboseLevel > 3)
    250242    G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
    251                                                                                                    
     243
    252244  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
    253245 
    254   if (photonEnergy0 <= LowEnergyLimit())
    255   {
     246  if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
     247    {
    256248      fParticleChange->ProposeTrackStatus(fStopAndKill);
    257249      fParticleChange->SetProposedKineticEnergy(0.);
    258250      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
    259251      return ;
    260   }
     252    }
    261253
    262254  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
    … …  
    534526  Xhigh=1e06;
    535527  if (samplingFunction_x)
    536     delete samplingFunction_x;
     528    {
     529      delete samplingFunction_x;
     530      samplingFunction_x = 0;
     531    }
    537532  if (samplingFunction_xNoLog)
    538     delete samplingFunction_xNoLog;
    539  
     533    {
     534      delete samplingFunction_xNoLog;
     535      samplingFunction_xNoLog = 0;
     536    }
     537
    540538  samplingFunction_x = new G4DataVector();
    541539  samplingFunction_xNoLog = new G4DataVector();

  • trunk/source/processes/electromagnetic/lowenergy/src/G4ShellVacancy.cc

    r819 r1055  
    5050
    5151{
     52/*
    5253  G4int size = xsis.size();
    5354  for (G4int k =0; k<size; k++)
    … …  
    5657      xsis[k] = 0;
    5758    }
     59*/
    5860}
    5961
    60 void G4ShellVacancy::AddXsiTable(G4VEMDataSet* set)
    61 
     62void G4ShellVacancy::AddXsiTable(G4VEMDataSet* p)
    6263{
    63   xsis.push_back(set);
     64  xsis.push_back(p);
    6465}
    6566

  • trunk/source/processes/electromagnetic/lowenergy/src/G4eIonisationCrossSectionHandler.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eIonisationCrossSectionHandler.cc,v 1.11 2006/06/29 19:42:00 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4eIonisationCrossSectionHandler.cc,v 1.12 2009/01/29 08:13:34 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4242// 19 Jul 2002   VI          Create composite data set for material
    4343// 21 Jan 2003  V.Ivanchenko Cut per region
     44// 28 Jan 2009  L.Pandola    Added public method to make a easier migration of
     45//                           G4LowEnergyIonisation to G4LivermoreIonisationModel
    4446//
    4547// -------------------------------------------------------------------
    … …  
    98100    G4int nElements = material->GetNumberOfElements();
    99101
    100     if(verbose > 0) {
    101       G4cout << "eIonisation CS for " << m << "th material "
    102              << material->GetName()
    103              << "  eEl= " << nElements << G4endl;
    104     }
     102    if(verbose > 0)
     103      {
     104        G4cout << "eIonisation CS for " << m << "th material "
     105               << material->GetName()
     106               << "  eEl= " << nElements << G4endl;
     107      }
    105108
    106109    G4double tcut  = (*energyCuts)[m];
    … …  
    129132            value += cross * p * density;
    130133
    131             if(verbose>0 && m == 0 && e>=1. && e<=0.) {
     134            if(verbose>0 && m == 0 && e>=1. && e<=0.)
     135            {
    132136              G4cout << "G4eIonCrossSH: e(MeV)= " << e/MeV
    133137                     << " n= " << n
    … …  
    139143                     << " Z= " << Z
    140144                     << G4endl;
    141             }
     145              }
    142146
    143147          }
    … …  
    155159}
    156160
    157 
     161G4double G4eIonisationCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
     162                                                                                   G4double cutEnergy,
     163                                                                                   G4int Z)
     164{
     165  G4int nShells = NumberOfComponents(Z);
     166  G4double value = 0.;
     167  if(energy > cutEnergy)
     168    {
     169      for (G4int n=0; n<nShells; n++) {
     170        G4double cross = FindValue(Z, energy, n);
     171        G4double p = theParam->Probability(Z, cutEnergy, energy, energy, n);
     172        value += cross * p;
     173      }
     174    }
     175  return value;
     176}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4eIonisationParameters.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4eIonisationParameters.cc,v 1.23 2006/06/29 19:42:02 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4eIonisationParameters.cc,v 1.24 2009/03/23 09:13:44 pandola Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
    … …  
    4141//                          chenged to "ion-..."
    4242// 17.02.04 V.Ivanchenko    Increase buffer size
     43// 23.03.09 L.Pandola       Updated warning message
    4344//
    4445// -------------------------------------------------------------------
    … …  
    199200    if (! (lsdp->is_open()) ) {
    200201      G4String excep = G4String("G4IonisationParameters - data file: ")
    201       + name + G4String(" not found. The version 1.# of G4LEDATA should be used");
     202      + name + G4String(" not found. Please check and/or update G4LEDATA");
    202203      G4Exception(excep);
    203204    }

  • trunk/source/processes/electromagnetic/lowenergy/src/G4eIonisationSpectrum.cc

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eIonisationSpectrum.cc,v 1.25 2006/06/29 19:42:04 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4eIonisationSpectrum.cc,v 1.26 2009/03/23 09:13:44 pandola Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646// 30.05.2002 VI            Update to 24-parameters data
    4747// 11.07.2002 VI            Fix in integration over spectrum
     48// 23.03.2009 LP            Added protection against division by zero (for
     49//                          faulty database files), Bug Report 1042
    4850//
    4951// -------------------------------------------------------------------
    … …  
    125127  p.push_back((2.0*g - 1.0)/(g*g));
    126128 
    127   p[iMax-1] = Function(p[3], p);
     129  //Add protection against division by zero: actually in Function()
     130  //parameter p[3] appears in the denominator. Bug report 1042
     131  if (p[3] > 0)
     132    p[iMax-1] = Function(p[3], p);
     133  else
     134    {
     135      G4cout << "WARNING: G4eIonisationSpectrum::Probability "
     136             << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
     137             << Z << ". Please check and/or update it " << G4endl;     
     138    }
    128139
    129140  if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
    … …  
    210221  p.push_back((2.0*g - 1.0)/(g*g));
    211222
    212   p[iMax-1] = Function(p[3], p);
     223 
     224  //Add protection against division by zero: actually in Function()
     225  //parameter p[3] appears in the denominator. Bug report 1042
     226  if (p[3] > 0)
     227    p[iMax-1] = Function(p[3], p);
     228  else
     229    {
     230      G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
     231             << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
     232             << Z << ". Please check and/or update it " << G4endl;     
     233    }
    213234   
    214235  G4double val = AverageValue(x1, x2, p);
    … …  
    288309  p.push_back((2.0*g - 1.0)/(g*g));
    289310
    290   p[iMax-1] = Function(p[3], p);
     311 
     312  //Add protection against division by zero: actually in Function()
     313  //parameter p[3] appears in the denominator. Bug report 1042
     314  if (p[3] > 0)
     315    p[iMax-1] = Function(p[3], p);
     316  else
     317    {
     318      G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
     319             << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
     320             << Z << ". Please check and/or update it " << G4endl;     
     321    }
    291322
    292323  G4double aria1 = 0.0;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4hLowEnergyLoss.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4hLowEnergyLoss.cc,v 1.27 2008/06/20 19:54:03 sincerti Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4hLowEnergyLoss.cc,v 1.28 2009/02/20 10:49:54 sincerti Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// -----------------------------------------------------------
    … …  
    145145G4double G4hLowEnergyLoss::HighestKineticEnergy= 100.*GeV;
    146146G4int    G4hLowEnergyLoss::TotBin = 360;
    147 G4double G4hLowEnergyLoss::RTable,G4hLowEnergyLoss::LOGRTable;
     147G4double G4hLowEnergyLoss::RTable =1.1;
     148G4double G4hLowEnergyLoss::LOGRTable = 1.1;
    148149
    149150//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    772773    theRangeTable = theRangepbarTable ;
    773774  }
    774  
    775775  G4double R2 = RTable*RTable ;
    776776  G4double R1 = RTable+1.;
    … …  
    10621062{
    10631063  G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
    1064   G4double Tbin = LowestKineticEnergy/RTable ;
     1064  G4double Tbin = 0;
     1065  if (RTable !=0.) Tbin = LowestKineticEnergy/RTable ;
    10651066  G4double rangebin = 0.0 ;
    10661067  G4int binnumber = -1 ;
Note: See TracChangeset for help on using the changeset viewer.