Ignore:
Timestamp:
Jan 8, 2010, 3:02:48 PM (16 years ago)
Author:
garnier
Message:

update to geant4.9.3

Location:
trunk/examples/extended/electromagnetic/TestEm7/src
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/extended/electromagnetic/TestEm7/src/DetectorConstruction.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: DetectorConstruction.cc,v 1.8 2007/01/11 15:41:46 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: DetectorConstruction.cc,v 1.10 2008/04/21 13:13:30 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    4545
    4646#include "G4NistManager.hh"
    47 
    4847#include "G4UnitsTable.hh"
     48
     49#include "G4FieldManager.hh"
     50#include "G4TransportationManager.hh"
     51#include "G4RunManager.hh"
    4952
    5053//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    117120  Air->AddElement(N, fractionmass=0.7);
    118121  Air->AddElement(O, fractionmass=0.3);
     122
     123  density     = 1.e-5*g/cm3;
     124  pressure    = 2.e-2*bar;
     125  temperature = STP_Temperature;  // From PhysicalConstants.h .
     126  G4Material* vac = new G4Material( "TechVacuum", density, 1,
     127                           kStateGas, temperature, pressure );
     128  vac->AddMaterial( Air, 1. );
    119129
    120130  density     = universe_mean_density;    //from PhysicalConstants.h
    … …  
    235245{
    236246  absorSizeX = value; worldSizeX = 1.2*absorSizeX;
     247  G4RunManager::GetRunManager()->GeometryHasBeenModified();
    237248}
    238249 
    … …  
    243254  absorSizeYZ = value;
    244255  worldSizeYZ = 1.2*absorSizeYZ;
     256  G4RunManager::GetRunManager()->GeometryHasBeenModified();
    245257} 
    246258
    … …  
    252264  G4Material* pttoMaterial =
    253265    G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
    254   if (pttoMaterial) absorMaterial = pttoMaterial;
    255 }
    256 
    257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    258 
    259 #include "G4FieldManager.hh"
    260 #include "G4TransportationManager.hh"
     266  if (pttoMaterial) {
     267    absorMaterial = pttoMaterial;
     268    if(lAbsor) {
     269      lAbsor->SetMaterial(absorMaterial);
     270      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
     271    }
     272  }
     273}
     274
     275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    261276
    262277void DetectorConstruction::SetMagField(G4double fieldValue)
    … …  
    285300{
    286301  tallySize = value;
     302  G4RunManager::GetRunManager()->GeometryHasBeenModified();
    287303} 
    288304
    … …  
    294310  G4Material* pttoMaterial =
    295311    G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
    296   if (pttoMaterial) tallyMaterial = pttoMaterial;
     312  if (pttoMaterial) {
     313    tallyMaterial = pttoMaterial;
     314    if(lTally) {
     315      lTally->SetMaterial(tallyMaterial);
     316      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
     317    }
     318  }
    297319}
    298320
    … …  
    305327    tallyNumber++;
    306328  }
     329  G4RunManager::GetRunManager()->GeometryHasBeenModified();
    307330} 
    308331
    309332//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    310 
    311 #include "G4RunManager.hh"
    312333 
    313334void DetectorConstruction::UpdateGeometry()
    314335{
    315 G4RunManager::GetRunManager()->DefineWorldVolume(ConstructVolumes());
    316 }
    317 
    318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     336  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
     337  G4RunManager::GetRunManager()->DefineWorldVolume(ConstructVolumes());
     338}
     339
     340//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/DetectorMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: DetectorMessenger.cc,v 1.3 2006/06/29 16:58:13 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/EventAction.cc

    r807 r1230  
    2525//
    2626// $Id: EventAction.cc,v 1.4 2006/06/29 16:58:15 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/EventActionMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: EventActionMessenger.cc,v 1.4 2006/06/29 16:58:17 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/G4ScreenedNuclearRecoil.cc

    r807 r1230  
    2525//
    2626//
    27 // $Id: G4ScreenedNuclearRecoil.cc,v 1.5 2008/01/14 12:11:39 vnivanch Exp $
    28 // GEANT4 tag $Name:  $
     27// G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp
     28// GEANT4 tag
    2929//
    3030//
    … …  
    8484#include "G4ScreenedNuclearRecoil.hh"
    8585
     86const char* G4ScreenedCoulombCrossSectionInfo::CVSFileVers() { return
     87        "G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp GEANT4 tag ";
     88}
     89
    8690#include "G4ParticleTypes.hh"
    8791#include "G4ParticleTable.hh"
    8892#include "G4VParticleChange.hh"
    89 #include "G4ParticleChange.hh"
     93#include "G4ParticleChangeForLoss.hh"
    9094#include "G4DataVector.hh"
    9195#include "G4Track.hh"
    … …  
    99103#include "G4IsotopeVector.hh"
    100104
     105#include "G4EmProcessSubType.hh"
     106
    101107#include "G4RangeTest.hh"
    102108#include "G4ParticleDefinition.hh"
    … …  
    104110#include "G4ProcessManager.hh"
    105111#include "G4StableIsotopes.hh"
     112#include "G4LindhardPartition.hh"
    106113
    107114#include "Randomize.hh"
    … …  
    112119#include <iomanip>
    113120
     121#include "c2_factory.hh"
     122static c2_factory<G4double> c2; // this makes a lot of notation shorter
     123typedef c2_ptr<G4double> c2p;
     124
    114125G4ScreenedCoulombCrossSection::~G4ScreenedCoulombCrossSection()
    115126{
    116         ScreeningMap::iterator tables=screeningData.begin();
    117         for (;tables != screeningData.end(); tables++) {
    118                 delete (*tables).second.EMphiData;
    119         }
    120127        screeningData.clear();
    121        
    122         std::map<G4int, c2_function<G4double> *>::iterator mfpit=MFPTables.begin();
    123         for (;mfpit != MFPTables.end(); mfpit++) {
    124                 delete (*mfpit).second;
    125         }
    126         MFPTables.clear();
    127        
     128        MFPTables.clear();     
    128129}
    129130
    … …  
    248249                        element=elementVector[kel];
    249250                        G4int Z=(G4int)std::floor(element->GetZ()+0.5);
    250                         c2_function<G4double> &ifunc=*sigmaMap[Z];
     251                        const G4_c2_function &ifunc=sigmaMap[Z];
    251252                        if(!kel || ifunc.xmin() > emin) emin=ifunc.xmin();
    252253                        if(!kel || ifunc.xmax() < emax) emax=ifunc.xmax();                     
    … …  
    268269                        element=elementVector[kel];
    269270                        G4int Z=(G4int)std::floor(element->GetZ()+0.5);
    270                         c2_function<G4double> &sigma=*sigmaMap[Z];
     271                        const G4_c2_function &sigma=sigmaMap[Z];
    271272                        G4double ndens = atomDensities[kel]; // compute atom fraction for this element in this material
    272273                                               
    … …  
    281282                }
    282283                // and make a new interpolating function out of the sum
    283                 MFPTables[matidx] = static_cast<c2_function<G4double> *>(new log_log_interpolating_function<G4double>(
    284                                 evals, mfpvals));
     284                MFPTables[matidx] = c2.log_log_interpolating_function().load(evals, mfpvals,true,0,true,0);
    285285    }
    286        
    287 #ifdef DEBUG   
    288         for (G4int matidx=0; matidx < nMaterials; matidx++) {
    289                 const G4Material* material= (*materialTable)[matidx];
    290                 G4cout << "***** MFP (1MeV) ***** " << material->GetName() << "  " << (*MFPTables[matidx])(1.0) << G4endl;
    291         }
    292 #endif
    293        
    294286}
    295287
    … …  
    304296        recoilCutoff(RecoilCutoff), physicsCutoff(PhysicsCutoff),
    305297        hardeningFraction(0.0), hardeningFactor(1.0),
    306         externalCrossSectionConstructor(0)
    307 {
     298        externalCrossSectionConstructor(0),
     299        NIELPartitionFunction(new G4LindhardRobinsonPartition)
     300{
     301                // for now, point to class instance of this. Doing it by creating a new one fails
     302                // to correctly update NIEL
     303                // not even this is needed... done in G4VProcess().
     304                // pParticleChange=&aParticleChange;
     305                processMaxEnergy=50000.0*MeV;
    308306                highEnergyLimit=100.0*MeV;
    309307                lowEnergyLimit=physicsCutoff;
    … …  
    313311                AddStage(new G4ScreenedCoulombClassicalKinematics);
    314312                AddStage(new G4SingleScatter);
     313                SetProcessSubType(fCoulombScattering);
    315314}
    316315
    317316void G4ScreenedNuclearRecoil::ResetTables()
    318317{
    319         std::map<G4int, c2_function<G4double>*>::iterator xh=meanFreePathTables.begin();
    320         for(;xh != meanFreePathTables.end(); xh++) {
    321                 delete (*xh).second;
    322         }
    323         meanFreePathTables.clear();
    324318       
    325319        std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xt=crossSectionHandlers.begin();
    … …  
    338332
    339333        collisionStages.clear();
     334}
     335
     336void G4ScreenedNuclearRecoil::SetNIELPartitionFunction(const G4VNIELPartition *part)
     337{
     338        if(NIELPartitionFunction) delete NIELPartitionFunction;
     339        NIELPartitionFunction=part;
     340}
     341
     342void G4ScreenedNuclearRecoil::DepositEnergy(G4int z1, G4double a1, const G4Material *material, G4double energy)
     343{
     344        if(!NIELPartitionFunction) {
     345                IonizingLoss+=energy;
     346        } else {
     347                G4double part=NIELPartitionFunction->PartitionNIEL(z1, a1, material, energy);
     348                IonizingLoss+=energy*(1-part);
     349                NIEL += energy*part;
     350        }
    340351}
    341352
    … …  
    363374G4double G4ScreenedNuclearRecoil::GetMeanFreePath(const G4Track& track,
    364375                                                  G4double,
    365                                                   G4ForceCondition*)
     376                                                  G4ForceCondition* cond)
    366377{
    367378        const G4DynamicParticle* incoming = track.GetDynamicParticle();
    … …  
    370381       
    371382        G4double meanFreePath;
    372        
    373         if (energy < lowEnergyLimit || energy < recoilCutoff) return 1.0*nanometer; /* stop slow particles! */
    374         else if (energy > highEnergyLimit*a1) energy=highEnergyLimit*a1; /* constant MFP at high energy */
     383        *cond=NotForced;
     384       
     385        if (energy < lowEnergyLimit || energy < recoilCutoff*a1) {
     386                *cond=Forced;
     387                return 1.0*nm; /* catch and stop slow particles to collect their NIEL! */
     388        } else if (energy > processMaxEnergy*a1) {
     389                return DBL_MAX; // infinite mean free path
     390        } else if (energy > highEnergyLimit*a1) energy=highEnergyLimit*a1; /* constant MFP at high energy */
    375391
    376392        G4double fz1=incoming->GetDefinition()->GetPDGCharge();
    … …  
    390406        size_t materialIndex = materialCouple->GetMaterial()->GetIndex();
    391407
    392         c2_function<G4double> &mfp=*(*xs)[materialIndex];
     408        const G4_c2_function &mfp=*(*xs)[materialIndex];
    393409       
    394410        // make absolutely certain we don't get an out-of-range energy
    … …  
    403419{
    404420        validCollision=1;
    405         aParticleChange.Initialize(aTrack);
     421        pParticleChange->Initialize(aTrack);
    406422        NIEL=0.0; // default is no NIEL deposited
     423        IonizingLoss=0.0;
    407424       
    408425        // do universal setup
    … …  
    413430        G4double fz1=baseParticle->GetPDGCharge()/eplus;
    414431        G4int z1=(G4int)(fz1+0.5);
     432        G4double a1=baseParticle->GetPDGMass()/amu_c2;
    415433        G4double incidentEnergy = incidentParticle->GetKineticEnergy();
    416434               
    … …  
    418436        const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
    419437       
    420         if(incidentEnergy < GetRecoilCutoff()) { // check energy sanity on entry
    421                 if(!baseParticle->GetProcessManager()->
    422                    GetAtRestProcessVector()->size())
    423                         aParticleChange.ProposeTrackStatus(fStopAndKill);
    424                 else
    425                         aParticleChange.ProposeTrackStatus(fStopButAlive);
    426 
    427                 AddToNIEL(incidentEnergy);
    428                 aParticleChange.ProposeEnergy(0.0);
     438        const G4Material* mat = couple->GetMaterial();
     439
     440        G4double P=0.0; // the impact parameter of this collision
     441
     442        if(incidentEnergy < GetRecoilCutoff()*a1) { // check energy sanity on entry
     443                DepositEnergy(z1, baseParticle->GetPDGMass()/amu_c2, mat, incidentEnergy);
     444                GetParticleChange().ProposeEnergy(0.0);
    429445                // stop the particle and bail out
    430446                validCollision=0;
    431         }       
    432        
    433         const G4Material* mat = couple->GetMaterial();
    434         G4double numberDensity=mat->GetTotNbOfAtomsPerVolume();
    435         G4double lattice=0.5/std::pow(numberDensity,1.0/3.0); // typical lattice half-spacing
    436         G4double length=GetCurrentInteractionLength();
    437         G4double sigopi=1.0/(CLHEP::pi*numberDensity*length);  // this is sigma0/pi
    438        
    439         // compute the impact parameter very early, so if is rejected as too far away, little effort is wasted
    440         // this is the TRIM method for determining an impact parameter based on the flight path
    441         // this gives a cumulative distribution of N(P)= 1-exp(-pi P^2 n l)
    442         // which says the probability of NOT hitting a disk of area sigma= pi P^2 =exp(-sigma N l)
    443         // which may be reasonable
    444         G4double P;
    445         if(sigopi < lattice*lattice) {
    446                 // normal long-flight approximation
    447                 P = std::sqrt(-std::log(G4UniformRand()) *sigopi);
    448447        } else {
    449                 // short-flight limit
    450                 P = std::sqrt(G4UniformRand())*lattice;
    451         }
    452                
    453         G4double fraction=GetHardeningFraction();
    454         if(fraction && G4UniformRand() < fraction) {
    455                 // pick out some events, and increase the central cross section
    456                 // by reducing the impact parameter
    457                 P /= std::sqrt(GetHardeningFactor());
    458         }
    459        
    460        
    461         // check if we are far enough away that the energy transfer must be below cutoff,
    462         // and leave everything alone if so, saving a lot of time.
    463         if(P*P > sigopi) {
    464                 if(GetVerboseLevel() > 1)
    465                         printf("ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n",
    466                                    length/angstrom, P/angstrom,std::sqrt(sigopi)/angstrom);
    467                 // no collision, don't follow up with anything
    468                 validCollision=0;
     448               
     449                G4double numberDensity=mat->GetTotNbOfAtomsPerVolume();
     450                G4double lattice=0.5/std::pow(numberDensity,1.0/3.0); // typical lattice half-spacing
     451                G4double length=GetCurrentInteractionLength();
     452                G4double sigopi=1.0/(CLHEP::pi*numberDensity*length);  // this is sigma0/pi
     453               
     454                // compute the impact parameter very early, so if is rejected as too far away, little effort is wasted
     455                // this is the TRIM method for determining an impact parameter based on the flight path
     456                // this gives a cumulative distribution of N(P)= 1-exp(-pi P^2 n l)
     457                // which says the probability of NOT hitting a disk of area sigma= pi P^2 =exp(-sigma N l)
     458                // which may be reasonable
     459                if(sigopi < lattice*lattice) {
     460                        // normal long-flight approximation
     461                        P = std::sqrt(-std::log(G4UniformRand()) *sigopi);
     462                } else {
     463                        // short-flight limit
     464                        P = std::sqrt(G4UniformRand())*lattice;
     465                }
     466                       
     467                G4double fraction=GetHardeningFraction();
     468                if(fraction && G4UniformRand() < fraction) {
     469                        // pick out some events, and increase the central cross section
     470                        // by reducing the impact parameter
     471                        P /= std::sqrt(GetHardeningFactor());
     472                }
     473               
     474               
     475                // check if we are far enough away that the energy transfer must be below cutoff,
     476                // and leave everything alone if so, saving a lot of time.
     477                if(P*P > sigopi) {
     478                        if(GetVerboseLevel() > 1)
     479                                printf("ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n",
     480                                           length/angstrom, P/angstrom,std::sqrt(sigopi)/angstrom);
     481                        // no collision, don't follow up with anything
     482                        validCollision=0;
     483                }
    469484        }
    470485       
    471486        // find out what we hit, and record it in our kinematics block.
     487        kinematics.targetMaterial=mat;
     488        kinematics.a1=a1;
     489       
    472490        if(validCollision) {
    473491                G4ScreenedCoulombCrossSection   *xsect=GetCrossSectionHandlers()[z1];
    … …  
    477495                kinematics.recoilIon=recoilIon;
    478496                kinematics.impactParameter=P;
    479                 kinematics.a1=baseParticle->GetPDGMass()/amu_c2;
    480497                kinematics.a2=recoilIon->GetPDGMass()/amu_c2;
    481498        } else {
    482499                kinematics.recoilIon=0;
    483500                kinematics.impactParameter=0;
    484                 kinematics.a1=baseParticle->GetPDGMass()/amu_c2;
    485501                kinematics.a2=0;
    486502        }
    … …  
    492508       
    493509        if(registerDepositedEnergy) {
    494           aParticleChange.ProposeLocalEnergyDeposit(NIEL);
    495           aParticleChange.ProposeNonIonizingEnergyDeposit(NIEL);
    496         }
     510                pParticleChange->ProposeLocalEnergyDeposit(IonizingLoss+NIEL);
     511                pParticleChange->ProposeNonIonizingEnergyDeposit(NIEL);
     512                //MHM G4cout << "depositing energy, total = " << IonizingLoss+NIEL << " NIEL = " << NIEL << G4endl;
     513        }
     514
    497515        return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
     516}
     517
     518G4ScreenedCoulombClassicalKinematics::G4ScreenedCoulombClassicalKinematics() :
     519// instantiate all the needed functions statically, so no allocation is done at run time
     520// we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
     521// or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
     522// note that only the last of these gets deleted, since it owns the rest
     523phifunc(c2.const_plugin_function()),
     524xovereps(c2.linear(0., 0., 0.)), // will fill this in with the right slope at run time
     525diff(c2.quadratic(0., 0., 0., 1.)-xovereps*phifunc)
     526{
    498527}
    499528
    … …  
    520549        }
    521550       
    522         c2_function<G4double> &phiData=*(screen->EMphiData);
    523         // instantiate all the needed functions statically, so no allocation is done at run time
    524551        // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
    525552        // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
    526         static c2_plugin_function<G4double> phifunc;
    527         static c2_quadratic<G4double> xsq(0., 0., 0., 1.); //  x^2
    528         static c2_linear<G4double> xovereps(0., 0., 0.); // will fill this in with the right slope at run time
    529         static c2_function<G4double> &xphi=xovereps*phifunc;
    530         static c2_function<G4double> &diff=xsq-xphi;
    531        
    532553        xovereps.reset(0., 0.0, au/eps); // slope of x*au/eps term
    533         phifunc.set_function(phiData); // install interpolating table
    534        
     554        phifunc.set_function(&(screen->EMphiData.get())); // install interpolating table
    535555        G4double xx1, phip, phip2;
    536         G4int root_error;
    537        
    538         xx1=diff.find_root(phiData.xmin(), std::min(10*xx0*au,phiData.xmax()),
    539                                            std::min(xx0*au, phiData.xmax()), beta*beta*au*au, &root_error, &phip, &phip2)/au;
     556        G4int root_error;       
     557        xx1=diff->find_root(phifunc.xmin(), std::min(10*xx0*au,phifunc.xmax()),
     558                                           std::min(xx0*au, phifunc.xmax()), beta*beta*au*au, &root_error, &phip, &phip2)/au;
    540559       
    541560        if(root_error) {
    542561                G4cout << "Screened Coulomb Root Finder Error" << G4endl;
    543562                G4cout << "au " << au << " A " << A << " a1 " << a1 << " xx1 " << xx1 << " eps " << eps << " beta " << beta << G4endl;
    544                 G4cout << " xmin " << phiData.xmin() << " xmax " << std::min(10*xx0*au,phiData.xmax()) ;
    545                 G4cout << " f(xmin) " << phifunc(phiData.xmin()) <<  " f(xmax) " << phifunc(std::min(10*xx0*au,phiData.xmax())) ;
    546                 G4cout << " xstart " << std::min(xx0*au, phiData.xmax()) <<  " target " <<  beta*beta*au*au ;
     563                G4cout << " xmin " << phifunc.xmin() << " xmax " << std::min(10*xx0*au,phifunc.xmax()) ;
     564                G4cout << " f(xmin) " << phifunc(phifunc.xmin()) <<  " f(xmax) " << phifunc(std::min(10*xx0*au,phifunc.xmax())) ;
     565                G4cout << " xstart " << std::min(xx0*au, phifunc.xmax()) <<  " target " <<  beta*beta*au*au ;
    547566                G4cout << G4endl;
    548567                throw c2_exception("Failed root find");
    549568        }
    550569       
    551         phifunc.unset_function(); // throws an exception if used without setting again
    552                                                           // phiprime is scaled by one factor of au because phi is evaluated at (xx0*au),
     570        // phiprime is scaled by one factor of au because phi is evaluated at (xx0*au),
    553571        G4double phiprime=phip*au;
    554572       
    … …  
    564582                G4double x, ff;
    565583                x=xx1/xvals[k];
    566                 ff=1.0/std::sqrt(1.0-phiData(x*au)/(x*eps)-beta*beta/(x*x));
     584                ff=1.0/std::sqrt(1.0-phifunc(x*au)/(x*eps)-beta*beta/(x*x));
    567585                alpha+=weights[k]*ff;
    568586        }
    569587       
     588        phifunc.unset_function(); // throws an exception if used without setting again
     589
    570590        G4double thetac1=CLHEP::pi*beta*alpha/xx1; // complement of CM scattering angle
    571591        G4double sintheta=std::sin(thetac1); //note sin(pi-theta)=sin(theta)
    … …  
    626646        kin.eRecoil=eRecoil;
    627647       
    628         if(incidentEnergy-eRecoil < master->GetRecoilCutoff()) {
    629                 if(!baseParticle->GetProcessManager()->
    630                    GetAtRestProcessVector()->size())
    631                         aParticleChange.ProposeTrackStatus(fStopAndKill);
    632                 else
    633                         aParticleChange.ProposeTrackStatus(fStopButAlive);
     648        if(incidentEnergy-eRecoil < master->GetRecoilCutoff()*a1) {
    634649                aParticleChange.ProposeEnergy(0.0);
    635                 master->AddToNIEL(incidentEnergy-eRecoil);
     650                master->DepositEnergy(int(screen->z1), a1, kin.targetMaterial, incidentEnergy-eRecoil);
    636651        }
    637652       
    638         if(master->GetEnableRecoils() && eRecoil > master->GetRecoilCutoff()) {
     653        if(master->GetEnableRecoils() && eRecoil > master->GetRecoilCutoff() * kin.a2) {
    639654                kin.recoilIon=recoilIon;               
    640655        } else {
    641656                kin.recoilIon=0; // this flags no recoil to be generated
    642                 master->AddToNIEL(eRecoil) ;
     657                master->DepositEnergy(Z, A, kin.targetMaterial, eRecoil) ;
    643658        }       
    644659}
    … …  
    706721#include <vector>
    707722
    708 static c2_function<G4double> &ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
     723G4_c2_function &ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
    709724{
    710725        static const size_t ncoef=4;
    … …  
    729744       
    730745        *auval=au;
    731         return *static_cast<c2_function<G4double> *>(new lin_log_interpolating_function<G4double>(r, phi, false, phiprime0));
    732 }
    733 
    734 static c2_function<G4double> &MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
     746        return c2.lin_log_interpolating_function().load(r, phi, false, phiprime0,true,0);
     747}
     748
     749G4_c2_function &MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
    735750{       
    736751        static const size_t ncoef=3;
    … …  
    755770       
    756771        *auval=au;
    757         return *static_cast<c2_function<G4double> *>(new lin_log_interpolating_function<G4double>(r, phi, false, phiprime0));
    758 }
    759 
    760 static c2_function<G4double> &LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
     772        return c2.lin_log_interpolating_function().load(r, phi, false, phiprime0,true,0);
     773}
     774
     775G4_c2_function &LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
    761776{       
    762777//from Loftager, Besenbacher, Jensen & Sorensen
    … …  
    781796       
    782797        *auval=au;
    783         return *static_cast<c2_function<G4double> *>(new lin_log_interpolating_function<G4double>(r, phi, false, logphiprime0*phi[0]));
     798        return c2.lin_log_interpolating_function().load(r, phi, false, logphiprime0*phi[0],true,0);
     799}
     800
     801G4_c2_function &LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
     802{       
     803// hybrid of LJ and ZBL, uses LJ if x < 0.25*auniv, ZBL if x > 1.5*auniv, and
     804/// connector in between.  These numbers are selected so the switchover
     805// is very near the point where the functions naturally cross.
     806        G4double auzbl, aulj;
     807       
     808        c2p zbl=ZBLScreening(z1, z2, npoints, rMax, &auzbl);
     809        c2p lj=LJScreening(z1, z2, npoints, rMax, &aulj);
     810
     811        G4double au=(auzbl+aulj)*0.5;
     812        lj->set_domain(lj->xmin(), 0.25*au);
     813        zbl->set_domain(1.5*au,zbl->xmax());
     814       
     815        c2p conn=c2.connector_function(lj->xmax(), lj, zbl->xmin(), zbl, true,0);
     816        c2_piecewise_function_p<G4double> &pw=c2.piecewise_function();
     817        c2p keepit(pw);
     818        pw.append_function(lj);
     819        pw.append_function(conn);
     820        pw.append_function(zbl);
     821       
     822        *auval=au;
     823        keepit.release_for_return();
     824        return pw;
    784825}
    785826
    … …  
    791832        AddScreeningFunction("lj", LJScreening);       
    792833        AddScreeningFunction("mol", MoliereScreening); 
     834        AddScreeningFunction("ljzbl", LJZBLScreening); 
    793835}
    794836
    … …  
    855897                       
    856898                        G4double au;
    857                         c2_function<G4double> &screen=sfunc(z1, Z, 200, 50.0*angstrom, &au); // generate the screening data
    858 
     899                        G4_c2_ptr screen=sfunc(z1, Z, 200, 50.0*angstrom, &au); // generate the screening data
    859900                        G4ScreeningTables st;
    860                         st.EMphiData=&screen; // this is our phi table
     901       
     902                        st.EMphiData=screen; //save our phi table
    861903                        st.z1=z1; st.m1=a1; st.z2=Z; st.m2=a2; st.emin=recoilCutoff;
    862904                        st.au=au;
    … …  
    869911                        //this rearranges to phi(x0)/(x0*eps) = 2*theta/pi - theta^2/pi^2
    870912                       
    871                         c2_linear<G4double> c2au(0.0, 0.0, au);
    872                         c2_composed_function<G4double> phiau(screen, c2au); // build phi(x*au) for dimensionless phi
    873                         c2_linear<G4double> c2eps(0.0, 0.0, 0.0); // will store an appropriate eps inside this in loop
    874                         c2_ratio<G4double> x0func(phiau, c2eps); // this will be phi(x)/(x*eps) when c2eps is correctly set
     913                        c2_linear_p<G4double> &c2eps=c2.linear(0.0, 0.0, 0.0); // will store an appropriate eps inside this in loop
     914                        G4_c2_ptr phiau=screen(c2.linear(0.0, 0.0, au));
     915                        G4_c2_ptr x0func(phiau/c2eps); // this will be phi(x)/(x*eps) when c2eps is correctly set
     916                        x0func->set_domain(1e-6*angstrom/au, 0.9999*screen->xmax()/au); // needed for inverse function
     917                        // use the c2_inverse_function interface for the root finder... it is more efficient for an ordered
     918                        // computation of values.
     919                        G4_c2_ptr x0_solution(c2.inverse_function(x0func));
    875920                       
    876921                        G4double m1c2=a1*amu_c2;
    … …  
    897942                                G4double q=theta/pi;
    898943                                // G4cout << ee << " " << m1c2 << " " << gamma << " " << eps << " " << theta << " " << q << G4endl;
    899                                 G4double x0= x0func.find_root(1e-6*angstrom/au, 0.9999*screen.xmax()/au, 1.0, 2*q-q*q);
     944                                // old way using root finder
     945                                // G4double x0= x0func->find_root(1e-6*angstrom/au, 0.9999*screen.xmax()/au, 1.0, 2*q-q*q);
     946                                // new way using c2_inverse_function which caches useful information so should be a bit faster
     947                                // since we are scanning this in strict order.
     948                                G4double x0=0;
     949                                try {
     950                                        x0=x0_solution(2*q-q*q);
     951                                } catch(c2_exception e) {
     952                                        G4Exception(
     953                                                G4String("G4ScreenedNuclearRecoil: failure in inverse solution to generate MFP Tables: ")+e.what()
     954                                        );
     955                                }
    900956                                G4double betasquared=x0*x0 - x0*phiau(x0)/eps; 
    901957                                G4double sigma=pi*betasquared*au*au;
    … …  
    903959                                data[idx]=sigma;
    904960                        }
    905                        
    906961                        screeningData[Z]=st;
    907                         sigmaMap[Z] = static_cast<c2_function<G4double> *>(new log_log_interpolating_function<G4double>(
    908                                 energies, data));
     962                        sigmaMap[Z] = c2.log_log_interpolating_function().load(energies, data, true,0,true,0);
    909963                }
    910964        }

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysListEmLivermore.cc

    r807 r1230  
    2626//
    2727// $Id: PhysListEmLivermore.cc,v 1.2 2006/12/11 20:13:53 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysListEmPenelope.cc

    r807 r1230  
    2626//
    2727// $Id: PhysListEmPenelope.cc,v 1.1 2006/11/22 18:56:21 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysListEmStandard.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: PhysListEmStandard.cc,v 1.11 2008/01/14 12:11:39 vnivanch Exp $
    28 // GEANT4 tag $Name:  $
     26// $Id: PhysListEmStandard.cc,v 1.19 2008/11/20 20:34:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2928//
    3029//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    3938#include "G4PhotoElectricEffect.hh"
    4039
    41 #include "G4MultipleScattering.hh"
     40#include "G4eMultipleScattering.hh"
    4241#include "G4hMultipleScattering.hh"
    4342
    … …  
    5150
    5251#include "G4hIonisation.hh"
     52#include "G4hBremsstrahlung.hh"
     53#include "G4hPairProduction.hh"
     54
    5355#include "G4ionIonisation.hh"
    5456
    5557#include "G4EmProcessOptions.hh"
     58#include "G4MscStepLimitType.hh"
    5659
    5760//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    8689    } else if (particleName == "e-") {
    8790      //electron
    88       pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
    89       pmanager->AddProcess(new G4eIonisation,        -1, 2,2);
    90       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
     91      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
     92      pmanager->AddProcess(new G4eIonisation,         -1, 2, 2);
     93      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
    9194           
    9295    } else if (particleName == "e+") {
    9396      //positron
    94       pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
    95       pmanager->AddProcess(new G4eIonisation,        -1, 2,2);
    96       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
    97       pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,4);
     97      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
     98      pmanager->AddProcess(new G4eIonisation,         -1, 2, 2);
     99      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
     100      pmanager->AddProcess(new G4eplusAnnihilation,    0,-1, 4);
    98101     
    99102    } else if( particleName == "mu+" ||
    100103               particleName == "mu-"    ) {
    101104      //muon 
    102       pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
    103       pmanager->AddProcess(new G4MuIonisation,       -1, 2,2);
    104       pmanager->AddProcess(new G4MuBremsstrahlung,   -1, 3,3);
    105       pmanager->AddProcess(new G4MuPairProduction,   -1, 4,4);       
     105      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     106      pmanager->AddProcess(new G4MuIonisation,        -1, 2, 2);
     107      pmanager->AddProcess(new G4MuBremsstrahlung,    -1, 3, 3);
     108      pmanager->AddProcess(new G4MuPairProduction,    -1, 4, 4);       
     109             
     110    } else if( particleName == "proton" ||
     111               particleName == "pi-" ||
     112               particleName == "pi+"    ) {
     113      //proton 
     114      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     115      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
     116      pmanager->AddProcess(new G4hBremsstrahlung,     -1, 3, 3);
     117      pmanager->AddProcess(new G4hPairProduction,     -1, 4, 4);       
    106118     
    107     } else if( particleName == "alpha" || particleName == "GenericIon" ) {
    108       pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
    109       pmanager->AddProcess(new G4ionIonisation,      -1, 2,2);
     119    } else if( particleName == "alpha" ||
     120               particleName == "He3" ||
     121               particleName == "GenericIon" ) {
     122      //Ions
     123      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     124      pmanager->AddProcess(new G4ionIonisation,       -1, 2, 2);   
    110125
    111126    } else if ((!particle->IsShortLived()) &&
    … …  
    113128               (particle->GetParticleName() != "chargedgeantino")) {
    114129      //all others charged particles except geantino
    115       pmanager->AddProcess(new G4hMultipleScattering,-1,1,1);
    116       pmanager->AddProcess(new G4hIonisation,        -1,2,2);
     130      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     131      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
    117132    }
    118133  }
    119   G4EmProcessOptions opt;
    120   opt.SetStepFunction(0.2, 10*um);
    121   //  opt.SetSkin(1.);
    122   opt.SetMinEnergy(0.1*keV);
    123   opt.SetMaxEnergy(100.*GeV);
    124   opt.SetDEDXBinning(360);
    125   opt.SetLambdaBinning(360);
    126   opt.SetLinearLossLimit(1.e-6);
    127134}
    128135

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysListEmStandardNR.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: PhysListEmStandardNR.cc,v 1.1 2008/01/14 12:11:39 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: PhysListEmStandardNR.cc,v 1.3 2008/05/09 08:30:59 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5656#include "G4CoulombScattering.hh"
    5757
     58#include "G4DummyModel.hh"
     59
    5860//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5961
    … …  
    7476
    7577  G4ScreenedNuclearRecoil* nucr = new G4ScreenedNuclearRecoil();
     78  G4double energyLimit = 100.*MeV;
     79  nucr->SetMaxEnergyForScattering(energyLimit);
    7680
    7781  theParticleIterator->reset();
    … …  
    109113             
    110114    } else if (particleName == "alpha" || particleName == "He3") {
     115      G4hMultipleScattering* msc = new G4hMultipleScattering();
     116      G4DummyModel* dm = new G4DummyModel();
     117      dm->SetLowEnergyLimit(0.0);
     118      dm->SetHighEnergyLimit(energyLimit);
     119      msc->AddEmModel(0, dm);
     120      pmanager->AddProcess(msc, -1, 1,1);
     121
    111122      G4ionIonisation* ion = new G4ionIonisation();
    112123      ion->ActivateNuclearStopping(false);
    113       pmanager->AddProcess(ion, -1, 1, 1);
     124      pmanager->AddProcess(ion, -1, 2, 2);
     125
    114126      pmanager->AddDiscreteProcess(nucr);     
    115127
    116128    } else if (particleName == "GenericIon" ) {
     129      G4hMultipleScattering* msc = new G4hMultipleScattering();
     130      G4DummyModel* dm = new G4DummyModel();
     131      dm->SetLowEnergyLimit(0.0);
     132      dm->SetHighEnergyLimit(energyLimit);
     133      msc->AddEmModel(0, dm);
     134      pmanager->AddProcess(msc, -1, 1,1);
     135
    117136      G4ionIonisation* ion = new G4ionIonisation();
    118137      ion->ActivateNuclearStopping(false);
    119138      ion->SetStepFunction(0.1, um);
    120       pmanager->AddProcess(ion, -1, 1, 1);
     139      pmanager->AddProcess(ion, -1, 2, 2);
     140
    121141      pmanager->AddDiscreteProcess(nucr);     
    122142
    … …  
    124144               particleName == "deuteron" ||
    125145               particleName == "triton") {
     146      G4hMultipleScattering* msc = new G4hMultipleScattering();
     147      G4DummyModel* dm = new G4DummyModel();
     148      dm->SetLowEnergyLimit(0.0);
     149      dm->SetHighEnergyLimit(energyLimit);
     150      msc->AddEmModel(0, dm);
     151      pmanager->AddProcess(msc, -1, 1,1);
     152
    126153      G4hIonisation* hion = new G4hIonisation();
    127154      hion->SetFluctModel(new G4IonFluctuations());
    128155      hion->SetStepFunction(0.1, 10.*um);
    129156      hion->ActivateNuclearStopping(false);
    130       pmanager->AddProcess(hion, -1, 1,1);
     157      pmanager->AddProcess(hion, -1, 2, 2);
     158
    131159      pmanager->AddDiscreteProcess(nucr);     
    132160     

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysListEmStandardSS.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: PhysListEmStandardSS.cc,v 1.5 2008/01/14 12:11:39 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: PhysListEmStandardSS.cc,v 1.8 2008/11/16 19:17:39 maire Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    3838#include "G4PhotoElectricEffect.hh"
    3939
    40 #include "G4MultipleScattering.hh"
     40#include "G4CoulombScattering.hh"
    4141
    4242#include "G4eIonisation.hh"
    … …  
    5050#include "G4hIonisation.hh"
    5151#include "G4ionIonisation.hh"
    52 #include "G4ionGasIonisation.hh"
    53 #include "G4IonFluctuations.hh"
    54 #include "G4CoulombScattering.hh"
     52
     53#include "G4EmProcessOptions.hh"
    5554
    5655//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    8584    } else if (particleName == "e-") {
    8685      //electron
    87       pmanager->AddProcess(new G4eIonisation,        -1, 1,1);
    88       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 2,2);
    8986      pmanager->AddDiscreteProcess(new G4CoulombScattering);     
     87      pmanager->AddProcess(new G4eIonisation,        -1, 1, 1);
     88      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 2, 2);
    9089           
    9190    } else if (particleName == "e+") {
    9291      //positron
    93       pmanager->AddProcess(new G4eIonisation,        -1, 1,1);
    94       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 2,2);
    95       pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,3);
    9692      pmanager->AddDiscreteProcess(new G4CoulombScattering);     
     93      pmanager->AddProcess(new G4eIonisation,        -1, 1, 1);
     94      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 2, 2);
     95      pmanager->AddProcess(new G4eplusAnnihilation,   0,-1, 3);
    9796           
    9897    } else if (particleName == "mu+" ||
    9998               particleName == "mu-"    ) {
    100       //muon 
    101       pmanager->AddProcess(new G4MuIonisation,       -1, 1,1);
    102       pmanager->AddProcess(new G4MuBremsstrahlung,   -1, 2,2);
    103       pmanager->AddProcess(new G4MuPairProduction,   -1, 3,3);
    104       pmanager->AddDiscreteProcess(new G4CoulombScattering);     
     99      //muon
     100      pmanager->AddDiscreteProcess(new G4CoulombScattering);       
     101      pmanager->AddProcess(new G4MuIonisation,       -1, 1, 1);
     102      pmanager->AddProcess(new G4MuBremsstrahlung,   -1, 2, 2);
     103      pmanager->AddProcess(new G4MuPairProduction,   -1, 3, 3);
    105104             
    106     } else if (particleName == "alpha" || particleName == "He3") {
    107       G4ionIonisation* ion = new G4ionIonisation();
    108       ion->SetStepFunction(0.1, um);
    109       ion->ActivateNuclearStopping(false);
    110       pmanager->AddProcess(ion,  -1, 1,1);
    111       pmanager->AddDiscreteProcess(new G4CoulombScattering);     
    112 
    113     } else if (particleName == "GenericIon" ) {
    114       G4ionGasIonisation* ion = new G4ionGasIonisation();
    115       ion->ActivateNuclearStopping(false);
    116       ion->SetStepFunction(0.1, um);
    117       pmanager->AddProcess(ion,  -1, 1,1);
    118       G4CoulombScattering* cs = new G4CoulombScattering();
    119       cs->SetBuildTableFlag(false);
    120       pmanager->AddDiscreteProcess(cs);
    121      
     105    } else if (particleName == "alpha" ||
     106               particleName == "He3" ||
     107               particleName == "GenericIon") {
     108      pmanager->AddDiscreteProcess(new G4CoulombScattering);           
     109      pmanager->AddProcess(new G4ionIonisation,      -1, 1, 1);
     110     
    122111    } else if ((!particle->IsShortLived()) &&
    123112               (particle->GetPDGCharge() != 0.0) &&
    124113               (particle->GetParticleName() != "chargedgeantino")) {
    125114      //all others charged particles except geantino
    126       G4hIonisation* hion = new G4hIonisation();
    127       hion->SetStepFunction(0.1, 10.*um);
    128       hion->ActivateNuclearStopping(false);
    129       pmanager->AddProcess(hion,  -1,1,1);
    130115      pmanager->AddDiscreteProcess(new G4CoulombScattering);     
     116      pmanager->AddProcess(new G4hIonisation,        -1, 1, 1);
    131117    }
    132118  }
     119 
     120  // Em options
     121  //
     122  // Main options and setting parameters are shown here.
     123  // Several of them have default values.
     124  //
     125  G4EmProcessOptions emOptions;
     126 
     127  //physics tables
     128  //
     129  emOptions.SetMinEnergy(100*eV);       //default   
     130  emOptions.SetMaxEnergy(100*TeV);      //default 
     131  emOptions.SetDEDXBinning(12*20);      //default=12*7 
     132  emOptions.SetLambdaBinning(12*20);    //default=12*7
     133  emOptions.SetSplineFlag(true);        //default
     134     
     135  //energy loss
     136  //
     137  emOptions.SetStepFunction(0.2,50*um); //default=(0.2, 1*mm)   
     138  emOptions.SetLinearLossLimit(1.e-2);  //default
     139   
     140  //ionization
     141  //
     142  emOptions.SetSubCutoff(false);        //default   
    133143}
    134144

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysicsList.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: PhysicsList.cc,v 1.28 2008/01/14 12:11:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: PhysicsList.cc,v 1.36 2008/11/21 12:53:13 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    3535#include "PhysListEmStandard.hh"
    3636#include "PhysListEmStandardSS.hh"
    37 #include "PhysListEmStandardIG.hh"
    3837#include "PhysListEmStandardNR.hh"
    3938#include "PhysListEmLivermore.hh"
    … …  
    4241#include "G4EmStandardPhysics_option1.hh"
    4342#include "G4EmStandardPhysics_option2.hh"
     43#include "G4EmStandardPhysics_option3.hh"
     44
     45#include "G4DecayPhysics.hh"
    4446
    4547#include "G4HadronElasticPhysics.hh"
    … …  
    5052#include "G4IonBinaryCascadePhysics.hh"
    5153
    52 #include "G4EmProcessOptions.hh"
    53 
    5454#include "G4LossTableManager.hh"
    5555#include "G4UnitsTable.hh"
    … …  
    5757#include "G4ProcessManager.hh"
    5858#include "G4Decay.hh"
     59
     60#include "StepMax.hh"
     61
     62#include "G4IonFluctuations.hh"
     63#include "G4IonParametrisedLossModel.hh"
    5964
    6065//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    7984
    8085  // EM physics
    81   emPhysicsList = new PhysListEmStandard(emName = "standard");
    82 
     86  emPhysicsList = new G4EmStandardPhysics(1);
     87  emName = G4String("emstandard");
     88
     89  // Deacy physics and all particles
     90  decPhysicsList = new G4DecayPhysics();
    8391}
    8492
    … …  
    8997  delete pMessenger;
    9098  delete emPhysicsList;
    91   for(size_t i=0; i<hadronPhys.size(); i++) delete hadronPhys[i];
    92 }
    93 
    94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    95 
    96 // Bosons
    97 #include "G4ChargedGeantino.hh"
    98 #include "G4Geantino.hh"
    99 #include "G4Gamma.hh"
    100 
    101 // leptons
    102 #include "G4MuonPlus.hh"
    103 #include "G4MuonMinus.hh"
    104 #include "G4NeutrinoMu.hh"
    105 #include "G4AntiNeutrinoMu.hh"
    106 
    107 #include "G4Electron.hh"
    108 #include "G4Positron.hh"
    109 #include "G4NeutrinoE.hh"
    110 #include "G4AntiNeutrinoE.hh"
    111 
    112 // Hadrons
    113 #include "G4MesonConstructor.hh"
    114 #include "G4BaryonConstructor.hh"
    115 #include "G4IonConstructor.hh"
     99  delete decPhysicsList;
     100  for(size_t i=0; i<hadronPhys.size(); i++) {delete hadronPhys[i];}
     101}
    116102
    117103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    119105void PhysicsList::ConstructParticle()
    120106{
    121 // pseudo-particles
    122   G4Geantino::GeantinoDefinition();
    123   G4ChargedGeantino::ChargedGeantinoDefinition();
    124  
    125 // gamma
    126   G4Gamma::GammaDefinition();
    127  
    128 // leptons
    129   G4Electron::ElectronDefinition();
    130   G4Positron::PositronDefinition();
    131   G4MuonPlus::MuonPlusDefinition();
    132   G4MuonMinus::MuonMinusDefinition();
    133 
    134   G4NeutrinoE::NeutrinoEDefinition();
    135   G4AntiNeutrinoE::AntiNeutrinoEDefinition();
    136   G4NeutrinoMu::NeutrinoMuDefinition();
    137   G4AntiNeutrinoMu::AntiNeutrinoMuDefinition(); 
    138 
    139 // mesons
    140   G4MesonConstructor mConstructor;
    141   mConstructor.ConstructParticle();
    142 
    143 // barions
    144   G4BaryonConstructor bConstructor;
    145   bConstructor.ConstructParticle();
    146 
    147 // ions
    148   G4IonConstructor iConstructor;
    149   iConstructor.ConstructParticle();
    150 }
    151 
    152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     107  decPhysicsList->ConstructParticle();
     108}
     109
     110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     111
    153112void PhysicsList::ConstructProcess()
    154113{
    … …  
    160119  //
    161120  emPhysicsList->ConstructProcess();
     121  em_config.AddModels();
     122
     123  // decay physics list
     124  //
     125  decPhysicsList->ConstructProcess();
    162126 
    163127  // hadronic physics lists
    164   for(size_t i=0; i<hadronPhys.size(); i++) hadronPhys[i]->ConstructProcess();
    165  
    166   // decay process
    167   //
    168   G4Decay* fDecayProcess = new G4Decay();
    169 
    170   theParticleIterator->reset();
    171   while( (*theParticleIterator)() ){
    172     G4ParticleDefinition* particle = theParticleIterator->value();
    173     G4ProcessManager* pmanager = particle->GetProcessManager();
    174 
    175     if (fDecayProcess->IsApplicable(*particle) && !particle->IsShortLived()) {
    176 
    177       pmanager ->AddProcess(fDecayProcess);
    178 
    179       // set ordering for PostStepDoIt and AtRestDoIt
    180       pmanager ->SetProcessOrdering(fDecayProcess, idxPostStep);
    181       pmanager ->SetProcessOrdering(fDecayProcess, idxAtRest);
    182 
    183     }
     128  for(size_t i=0; i<hadronPhys.size(); i++) {
     129    hadronPhys[i]->ConstructProcess();
    184130  }
    185131 
    186132  // step limitation (as a full process)
    187133  // 
    188   AddStepMax();
    189  
    190   G4EmProcessOptions opt;
    191   opt.SetDEDXBinning(480);
     134  AddStepMax(); 
    192135}
    193136
    … …  
    196139void PhysicsList::AddPhysicsList(const G4String& name)
    197140{
    198   if (verboseLevel>-1) {
     141  if (verboseLevel>1) {
    199142    G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
    200143  }
    … …  
    202145  if (name == emName) return;
    203146
    204   if (name == "standard") {
     147  if (name == "standard_local") {
    205148
    206149    emName = name;
    … …  
    212155    emName = name;
    213156    delete emPhysicsList;
    214     emPhysicsList = new G4EmStandardPhysics();
     157    emPhysicsList = new G4EmStandardPhysics(1);
    215158
    216159  } else if (name == "emstandard_opt1") {
    … …  
    225168    delete emPhysicsList;
    226169    emPhysicsList = new G4EmStandardPhysics_option2();
     170   
     171  } else if (name == "emstandard_opt3") {
     172
     173    emName = name;
     174    delete emPhysicsList;
     175    emPhysicsList = new G4EmStandardPhysics_option3();   
    227176
    228177  } else if (name == "standardSS") {
    … …  
    238187    emPhysicsList = new PhysListEmStandardNR(name);
    239188
    240   } else if (name == "standardIG") {
    241 
    242     emName = name;
    243     delete emPhysicsList;
    244     emPhysicsList = new PhysListEmStandardIG(name);
     189  } else if (name == "standardICRU73") {
     190
     191    emName = name;
     192    delete emPhysicsList;
     193    emPhysicsList = new PhysListEmStandard(name);
     194    em_config.SetExtraEmModel("GenericIon","ionIoni",
     195                              new G4IonParametrisedLossModel(),
     196                              "",0.0, 100.0*TeV,
     197                              new G4IonFluctuations());
     198    G4cout << "standardICRU73" << G4endl;
    245199
    246200  } else if (name == "livermore") {
    … …  
    288242//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    289243
    290 #include "StepMax.hh"
    291 
    292244void PhysicsList::AddStepMax()
    293245{
    … …  
    297249  theParticleIterator->reset();
    298250  while ((*theParticleIterator)()){
    299       G4ParticleDefinition* particle = theParticleIterator->value();
    300       G4ProcessManager* pmanager = particle->GetProcessManager();
    301 
    302       if (stepMaxProcess->IsApplicable(*particle) && pmanager)
    303         {
    304           pmanager ->AddDiscreteProcess(stepMaxProcess);
    305         }
     251    G4ParticleDefinition* particle = theParticleIterator->value();
     252    G4ProcessManager* pmanager = particle->GetProcessManager();
     253
     254    if (stepMaxProcess->IsApplicable(*particle) && pmanager)
     255      {
     256        pmanager ->AddDiscreteProcess(stepMaxProcess);
     257      }
    306258  }
    307259}

  • trunk/examples/extended/electromagnetic/TestEm7/src/PhysicsListMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: PhysicsListMessenger.cc,v 1.3 2006/06/29 16:58:37 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/PrimaryGeneratorAction.cc

    r807 r1230  
    2525//
    2626// $Id: PrimaryGeneratorAction.cc,v 1.2 2006/06/29 16:58:39 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/PrimaryGeneratorMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: PrimaryGeneratorMessenger.cc,v 1.3 2006/06/29 16:58:41 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/RunAction.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: RunAction.cc,v 1.21 2008/01/14 12:11:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: RunAction.cc,v 1.24 2008/08/22 18:30:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    4242
    4343#include "Randomize.hh"
    44 
    45 #ifdef G4ANALYSIS_USE
    46 #include "AIDA/AIDA.h"
    47 #endif
     44#include "Histo.hh"
    4845
    4946//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5148RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
    5249                     PrimaryGeneratorAction* kin)
    53 :detector(det), physics(phys), kinematic(kin), af(0), tree(0)
     50:detector(det), physics(phys), kinematic(kin)
    5451{
    5552  tallyEdep = new G4double[MaxTally];
    5653  binLength = offsetX = 0.;
    57   histo[0] = 0;
    58  
    59 #ifdef G4ANALYSIS_USE
    60   // Creating the analysis factory
    61   af = AIDA_createAnalysisFactory();
    62   if(!af) {
    63     G4cout << "RunAction::RunAction() :"
    64            << " problem creating the AIDA analysis factory."
    65            << G4endl;
    66   }       
    67 #endif 
     54  histo = new Histo();
     55  histo->setFileName("testem7");
     56  histo->add1D("1","Edep (MeV/mm) along absorber (mm)", 100, 0, 100);
     57  histo->add1D("2","Edep (MeV/mm) along absorber zoomed (mm)", 100, 0, 100);
     58  histo->add1D("3","Projectile range (mm)", 100, 0, 100);
    6859}
    6960
    … …  
    7364{
    7465  delete [] tallyEdep;
    75  
    76 #ifdef G4ANALYSIS_USE
    77   delete af; 
    78 #endif     
    79 }
    80 
    81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    82 
    83 void RunAction::bookHisto()
    84 {
    85   length  = detector->GetAbsorSizeX();
    86   G4double stepMax = physics->GetStepMaxProcess()->GetMaxStep();
    87   const G4int nbmin = 100;
    88   G4int nbBins = (int)(0.5 + length/stepMax);
    89   if (nbBins < nbmin) nbBins = nbmin;
    90   binLength = length/nbBins;
    91   offsetX   = 0.5*length;
    92  
    93 #ifdef G4ANALYSIS_USE
    94   if (!af) return;
    95 
    96   // Create a tree mapped to an hbook file.
    97   G4bool readOnly  = false;
    98   G4bool createNew = true;
    99   G4String options = "--noErrors uncompress";
    100   AIDA::ITreeFactory* tf  = af->createTreeFactory(); 
    101   tree = tf->create("testem7.hbook","hbook", readOnly, createNew, options);
    102   //tree = tf->create("testem7.root", "root",readOnly, createNew, options);
    103   //tree = tf->create("testem7.XML" , "XML" ,readOnly, createNew, options);
    104   delete tf;
    105   if (!tree) {
    106     G4cout << "RunAction::bookHisto()" << G4endl;
    107     return;
    108   }
    109 
    110   // Create a histogram factory, whose histograms will be handled by the tree
    111   AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
    112  
    113   // Create histogram
    114   histo[0] = hf->createHistogram1D("1","Edep (MeV/mm) along absorber (mm)",
    115              nbBins, 0, length/mm);
    116              
    117   delete hf;         
    118   G4cout << "\n----> Histogram Tree opened" << G4endl;
    119 #endif
    120 }
    121 
    122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    123 
    124 void RunAction::cleanHisto()
    125 {
    126 #ifdef G4ANALYSIS_USE
    127   tree->commit();       // Writing the histograms to the file
    128   tree->close();        // and closing the tree (and the file)
    129   delete tree;
    130   tree = 0;
    131  
    132   G4cout << "\n----> Histogram Tree saved" << G4endl; 
    133 #endif
     66  delete histo; 
    13467}
    13568
    … …  
    13871void RunAction::FillHisto(G4int ih, G4double x, G4double weight)
    13972{
    140 #ifdef G4ANALYSIS_USE
    141   if(histo[ih]) histo[ih]->fill(x, weight);
    142 #endif
     73  histo->fill(ih, x, weight);
    14374}
    14475
    … …  
    15687  //
    15788  nPrimarySteps = 0;
     89  nRange = 0;
    15890  projRange = projRange2 = 0.;
    15991  edeptot = eniel = 0.;
    16092  for (G4int j=0; j<MaxTally; j++) tallyEdep[j] = 0.;
    16193  kinematic->ResetEbeamCumul();
    162   bookHisto();     
     94
     95  // define "1" histogram binning
     96  length  = detector->GetAbsorSizeX();
     97  G4double stepMax = physics->GetStepMaxProcess()->GetMaxStep();
     98  const G4int nbmin = 100;
     99  G4int nbBins = (G4int)(0.5 + length/stepMax);
     100  if (nbBins < nbmin) nbBins = nbmin;
     101  binLength = length/nbBins;
     102  offsetX   = 0.5*length;
     103   
     104  // histogram "1" is defined by the length of the target
     105  // zoomed histograms are defined by UI command
     106  histo->setHisto1D(0, nbBins, 0, length, mm);
     107 
     108  histo->book();
    163109}
    164110
    … …  
    186132  //compute projected range and straggling
    187133  //
    188   projRange /= NbofEvents; projRange2 /= NbofEvents;
     134  if(nRange > 0) {
     135    projRange /= nRange;
     136    projRange2 /= nRange;
     137  }
    189138  G4double rms = projRange2 - projRange*projRange;       
    190139  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
    … …  
    227176  }
    228177
    229 #ifdef G4ANALYSIS_USE
    230178  // normalize histogram
    231179  G4double fac = (mm/MeV)/(NbofEvents *  binLength);
    232   histo[0]->scale(fac);
    233 #endif
    234  
     180  for (G4int j=0; j<3; j++) {histo->scale(j, fac);}
    235181 
    236182  // save and clean histo
    237   cleanHisto();
     183  histo->save();
    238184 
    239185  // show Rndm status

  • trunk/examples/extended/electromagnetic/TestEm7/src/StepMax.cc

    r807 r1230  
    2525//
    2626// $Id: StepMax.cc,v 1.5 2006/06/29 16:58:45 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/StepMaxMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: StepMaxMessenger.cc,v 1.3 2006/06/29 16:58:47 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/SteppingAction.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: SteppingAction.cc,v 1.12 2008/01/14 12:11:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: SteppingAction.cc,v 1.14 2008/08/22 18:30:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5151void SteppingAction::UserSteppingAction(const G4Step* aStep)
    5252{
    53  G4double edep = aStep->GetTotalEnergyDeposit();
    54  if (edep <= 0.) return;
     53  G4double edep = aStep->GetTotalEnergyDeposit();
     54  if (edep <= 0.) return;
    5555
    56  // G4cout << "edep= " << edep << "NIEL= " << aStep->GetNonIonizingEnergyDeposit()<<G4endl;
     56  // G4cout << "edep= " << edep << "NIEL= " << aStep->GetNonIonizingEnergyDeposit()<<G4endl;
    5757
    58  runAction->FillEdep(edep,aStep->GetNonIonizingEnergyDeposit());
     58  runAction->FillEdep(edep,aStep->GetNonIonizingEnergyDeposit());
    5959
    60  if(aStep->GetTrack()->GetTrackID() == 1) runAction->AddPrimaryStep(); 
     60  if(aStep->GetTrack()->GetTrackID() == 1) runAction->AddPrimaryStep(); 
    6161 
    62  //Bragg curve
    63  //     
    64  G4StepPoint* prePoint  = aStep->GetPreStepPoint();
    65  G4StepPoint* postPoint = aStep->GetPostStepPoint();
     62  //Bragg curve
     63  //   
     64  G4StepPoint* prePoint  = aStep->GetPreStepPoint();
     65  G4StepPoint* postPoint = aStep->GetPostStepPoint();
    6666   
    67  G4double x1 = prePoint->GetPosition().x(), x2 = postPoint->GetPosition().x(); 
    68  G4double x = runAction->GetOffsetX() + x1 + G4UniformRand()*(x2-x1);
    69  runAction->FillHisto(0, x/mm , edep);
     67  G4double x1 = prePoint->GetPosition().x(), x2 = postPoint->GetPosition().x(); 
     68  G4double x = runAction->GetOffsetX() + x1 + G4UniformRand()*(x2-x1);
     69  runAction->FillHisto(0, x/mm , edep);
     70  runAction->FillHisto(1, x/mm , edep);
    7071
    71  //fill tallies
    72  //
    73  G4TouchableHandle touchable = prePoint->GetTouchableHandle();
    74  G4LogicalVolume* lVolume = touchable->GetVolume()->GetLogicalVolume();
    75  if (lVolume == detector->GetLogicalTally())
    76      runAction->FillTallyEdep(touchable->GetCopyNumber(), edep);
     72  //fill tallies
     73  //
     74  G4TouchableHandle touchable = prePoint->GetTouchableHandle();
     75  G4LogicalVolume* lVolume = touchable->GetVolume()->GetLogicalVolume();
     76  if (lVolume == detector->GetLogicalTally())
     77    runAction->FillTallyEdep(touchable->GetCopyNumber(), edep);
    7778}
    7879

  • trunk/examples/extended/electromagnetic/TestEm7/src/SteppingVerbose.cc

    r807 r1230  
    2525//
    2626// $Id: SteppingVerbose.cc,v 1.3 2006/06/29 16:58:51 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm7/src/TrackingAction.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: TrackingAction.cc,v 1.3 2006/11/22 17:58:11 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: TrackingAction.cc,v 1.5 2008/08/22 18:30:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5050    G4double x = aTrack->GetPosition().x() + runAction->GetOffsetX();
    5151    if(x > runAction->GetLength()) x = runAction->GetLength();
    52     runAction->AddProjRange(x);
     52    //G4cout << " range= " << x << " x= " << aTrack->GetPosition().x()
     53    //     << " ofset= " << runAction->GetOffsetX() << G4endl;
     54    if(x > 0.0) runAction->AddProjRange(x);
     55    runAction->FillHisto(2, x/mm, 1.0);
    5356  } 
    5457}
Note: See TracChangeset for help on using the changeset viewer.