9.5算是升级好了，开始往9.6升级 - Githubissues

mirguest / MirguestIssueReport

5 stars 3 forks source link

9.5算是升级好了，开始往9.6升级 #86

Open mirguest opened 11 years ago

mirguest commented 11 years ago

发现，在9.6中，之前链接有问题的情况竟然没有了。

mirguest commented 11 years ago

出现的问题：

Added G4MuonNuclearProcess, meant to replace the old process G4MuNuclearInteraction by separating model and cross-section classes.

mirguest commented 11 years ago

9.6中需要研究这个了。

mirguest commented 11 years ago

看看，9.5中，G4MuNuclearInteraction是什么样的。顺便比较一下。

mirguest commented 11 years ago

https://savannah.cern.ch/bugs/?96913 Op Boundary的 SetModel 可以删除。

dengzy commented 11 years ago

DsPhysConsHadron.cc G4UHadronElasticProcess替换成 G4WHadronElasticProcess

dengzy commented 11 years ago

DetSim/PMTSim/src/dywPMTOpticalModel.cc中，修改了

    @@ -320,12 +290,12 @@ dywPMTOpticalModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
         {
           _photon_energy= energy;
           _wavelength= twopi*hbarc / energy;
    -      n_glass= _rindex_glass->Value( energy );
    +      n_glass= _rindex_glass->GetProperty( energy );
           _n1= n_glass; // just in case we exit before setting _n1
    -      _n2= _rindex_photocathode->Value( energy );
    -      _k2= _kindex_photocathode->Value( energy );
    +      _n2= _rindex_photocathode->GetProperty( energy );
    +      _k2= _kindex_photocathode->GetProperty( energy );
           _n3= 1.0;     // just in case we exit before setting _n3
    -      _efficiency= _efficiency_photocathode->Value( energy );
    +      _efficiency= _efficiency_photocathode->GetProperty( energy );
         }
       // initialize "whereAmI"
       if ( fastTrack.GetPrimaryTrack()->GetVolume() == _inner1_phys )
    @@ -403,7 +373,7 @@ dywPMTOpticalModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
           norm *= -1.0; // in principle, this is more efficient than norm= -norm;

         // set _thickness and _cos_theta1
    -    _thickness=  _thickness_photocathode->Value( pos.z() );
    +    _thickness=  _thickness_photocathode->GetProperty( pos.z() );
         _cos_theta1= dir * norm;

dengzy commented 11 years ago

DetSim/PhysiSim/include/DsG4Cerenkov.h

             G4double GetAverageNumberOfPhotons(const G4double charge,
                                                const G4double beta,
                                                const G4Material *aMaterial,
    -                                           G4MaterialPropertyVector* Rindex) const;
    +                                           const G4MaterialPropertyVector* Rindex) const;

dengzy commented 11 years ago

值得怀疑。

    diff --git a/DetSim/PhysiSim/include/DsG4GdNeutronHPCaptureFS.h b/DetSim/PhysiSim/include/DsG4GdNeutronHPCaptureFS.h
    index 3c7b196..1b31cd4 100644
    --- a/DetSim/PhysiSim/include/DsG4GdNeutronHPCaptureFS.h
    +++ b/DetSim/PhysiSim/include/DsG4GdNeutronHPCaptureFS.h
    @@ -35,7 +35,7 @@ public:
         ~DsG4GdNeutronHPCaptureFS();

         void   UpdateNucleus( const G4Fragment* , G4double );
    -    void Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType);
    +    void Init (G4double A, G4double Z, G4String & dirName, G4String & aFSType);

         G4HadFinalState * ApplyYourself(const G4HadProjectile & theTrack);
         G4NeutronHPFinalState * New() {

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/include/DsG4NNDCNeutronHPCaptureFS.h b/DetSim/PhysiSim/include/DsG4NNDCNeutronHPCaptureFS.h
    index 5d8ec26..f5ee2c4 100644
    --- a/DetSim/PhysiSim/include/DsG4NNDCNeutronHPCaptureFS.h
    +++ b/DetSim/PhysiSim/include/DsG4NNDCNeutronHPCaptureFS.h
    @@ -33,7 +33,7 @@ public:
         ~DsG4NNDCNeutronHPCaptureFS();

         void   UpdateNucleus( const G4Fragment* , G4double );
    -    void Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType);
    +    void Init (G4double A, G4double Z, G4String & dirName, G4String & aFSType);

         G4HadFinalState * ApplyYourself(const G4HadProjectile & theTrack);
         G4NeutronHPFinalState * New() {

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4Cerenkov.cc b/DetSim/PhysiSim/src/DsG4Cerenkov.cc
    index 0c97d23..896ffb5 100644
    --- a/DetSim/PhysiSim/src/DsG4Cerenkov.cc
    +++ b/DetSim/PhysiSim/src/DsG4Cerenkov.cc

    -        G4double Pmin = Rindex->GetMinLowEdgeEnergy();
    -        G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
    -        G4double dp = Pmax - Pmin;
    + G4double Pmin = Rindex->GetMinPhotonEnergy();
    + G4double Pmax = Rindex->GetMaxPhotonEnergy();
    + G4double dp = Pmax - Pmin;

    -        G4double nMax = Rindex->GetMaxValue();
    + G4double nMax = Rindex->GetMaxProperty();

dengzy commented 11 years ago

下面是迭代器变为for循环。注意！！！

    diff --git a/DetSim/PhysiSim/src/DsG4Cerenkov.cc b/DetSim/PhysiSim/src/DsG4Cerenkov.cc
    index 0c97d23..896ffb5 100644
    --- a/DetSim/PhysiSim/src/DsG4Cerenkov.cc
    +++ b/DetSim/PhysiSim/src/DsG4Cerenkov.cc
    @@ -459,97 +459,104 @@ DsG4Cerenkov::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
    +      if (theRefractionIndexVector) {
    +   
    +         // Retrieve the first refraction index in vector
    +         // of (photon energy, refraction index) pairs 

    -                      G4double currentRI = (*theRefractionIndexVector)[0];
    +         theRefractionIndexVector->ResetIterator();
    +         ++(*theRefractionIndexVector);  // advance to 1st entry 

    -                      if (currentRI > 1.0) {
    +         G4double currentRI = theRefractionIndexVector->
    +              GetProperty();

    -                         // Create first (photon energy, Cerenkov Integral)
    -                         // pair  
    +         if (currentRI > 1.0) {

    -                         G4double currentPM = theRefractionIndexVector->Energy(0);
    -                         G4double currentCAI = 0.0;

dengzy commented 11 years ago


    -                         for(G4int i = 1; i < theRefractionIndexVector->GetVectorLength(); ++i)
    -                         {
    -                                currentRI = (*theRefractionIndexVector)[i];    
    +      // loop over all (photon energy, refraction index)
    +      // pairs stored for this material  

    -                                currentPM = theRefractionIndexVector->Energy(i);
    +      while(++(*theRefractionIndexVector))
    +      {
    +       currentRI=theRefractionIndexVector->  
    +           GetProperty();

    -                                currentCAI = 0.5*(1.0/(prevRI*prevRI) +
    -                                                  1.0/(currentRI*currentRI));
    +       currentPM = theRefractionIndexVector->
    +           GetPhotonEnergy();

此处，i从1开始？对否？

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4OpBoundaryProcess.cc b/DetSim/PhysiSim/src/DsG4OpBoundaryProcess.cc
    index a3d1d8b..7dde5ba 100644
    --- a/DetSim/PhysiSim/src/DsG4OpBoundaryProcess.cc
    +++ b/DetSim/PhysiSim/src/DsG4OpBoundaryProcess.cc
             if (Rindex) {
    -                Rindex1 = Rindex->Value(thePhotonMomentum);
    -        }
    -        else {
    -                theStatus = NoRINDEX;
    -                aParticleChange.ProposeTrackStatus(fStopAndKill);
    -                return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
    -        }
    +   Rindex1 = Rindex->GetProperty(thePhotonMomentum);
    + }
    + else {
    +         theStatus = NoRINDEX;
    +   aParticleChange.ProposeTrackStatus(fStopAndKill);
    +   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
    + }

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    index 1f4e10d..e9140f2 100644
    --- a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    +++ b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    @@ -290,7 +290,7 @@ G4double DsG4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
                        aMaterialPropertyTable->GetProperty("RAYLEIGH");
                  if(AttenuationLengthVector){
                    AttenuationLength = AttenuationLengthVector ->
    -                                    Value(thePhotonEnergy);
    +                                    GetProperty(thePhotonEnergy);
                  }
                  else{
     //               G4cout << "No Rayleigh scattering length specified" << G4endl;

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    index 1f4e10d..e9140f2 100644
    --- a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    +++ b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    @@ -337,22 +337,24 @@ DsG4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aM
             G4double refraction_index;

             G4PhysicsOrderedFreeVector *RayleighScatteringLengths =
    -                                new G4PhysicsOrderedFreeVector();
    +       new G4PhysicsOrderedFreeVector();

             if (Rindex ) {

    -            for (G4int i=0; i < Rindex->GetVectorLength(); ++i) {
    +           Rindex->ResetIterator();

    -                e = (Rindex->Energy(i));
    +           while (++(*Rindex)) {

    -                refraction_index = (*Rindex)[i];
    +                e = (Rindex->GetPhotonEnergy());
    +
    +                refraction_index = Rindex->GetProperty();
                     refsq = refraction_index*refraction_index;
                     xlambda = h_Planck*c_light/e;

    -                if (verboseLevel>0) {
    -                        G4cout << Rindex->GetEnergy(i) << " MeV\t";
    -                        G4cout << xlambda << " mm\t";
    -                }
    +         if (verboseLevel>0) {
    +                 G4cout << Rindex->GetPhotonEnergy() << " MeV\t";
    +                 G4cout << xlambda << " mm\t";
    +   }

此处，为何又从零开始？

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    index 1f4e10d..e9140f2 100644
    --- a/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    +++ b/DetSim/PhysiSim/src/DsG4OpRayleigh.cc
    @@ -361,14 +363,14 @@ DsG4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aM

                     Dist = 1.0 / (c1*c2*c3*c4);

    -                if (verboseLevel>0) {
    -                        G4cout << Dist << " mm" << G4endl;
    -                }
    +         if (verboseLevel>0) {
    +                 G4cout << Dist << " mm" << G4endl;
    +   }
                     RayleighScatteringLengths->
    -                        InsertValues(Rindex->Energy(i), Dist);
    +     InsertValues(Rindex->GetPhotonEnergy(), Dist);
                }

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsG4Scintillation.cc b/DetSim/PhysiSim/src/DsG4Scintillation.cc
    index 446478b..6e88e7a 100644
    --- a/DetSim/PhysiSim/src/DsG4Scintillation.cc
    +++ b/DetSim/PhysiSim/src/DsG4Scintillation.cc
    @@ -265,7 +265,7 @@ DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
             aMaterialPropertiesTable->GetProperty("FASTCOMPONENT"); 
         const G4MaterialPropertyVector* Slow_Intensity =
             aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
    -    G4MaterialPropertyVector* Reemission_Prob =
    +    const G4MaterialPropertyVector* Reemission_Prob =
             aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
         if (verboseLevel > 0 ) {
           G4cout << " MaterialPropertyVectors: Fast_Intensity " << Fast_Intensity 
    @@ -296,7 +296,7 @@ DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
             if ( Reemission_Prob == 0)
                 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
             G4double p_reemission=
    -            Reemission_Prob->Value(aTrack.GetKineticEnergy());
    +            Reemission_Prob->GetProperty(aTrack.GetKineticEnergy());
             if (G4UniformRand() >= p_reemission)
                 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
             NumTracks= 1;
    @@ -322,22 +322,22 @@ DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)

             G4double ScintillationYield = 0;
             {// Yield.  Material must have this or we lack raisins dayetras
    -            G4MaterialPropertyVector* ptable =
    +            const G4MaterialPropertyVector* ptable =
                     aMaterialPropertiesTable->GetProperty("SCINTILLATIONYIELD");
                 if (!ptable) {
                     G4cout << "ConstProperty: failed to get SCINTILLATIONYIELD"
                            << G4endl;
                     return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
                 }
    -            ScintillationYield = ptable->Value(0);
    +            ScintillationYield = ptable->GetProperty(0);
             }

dengzy commented 11 years ago

    @@ -772,7 +772,11 @@ void DsG4Scintillation::BuildThePhysicsTable()
                     // Retrieve the first intensity point in vector
                     // of (photon energy, intensity) pairs

    -                G4double currentIN = (*theFastLightVector)[0];
    +                theFastLightVector->ResetIterator();
    +                ++(*theFastLightVector);        // advance to 1st entry 
    +
    +                G4double currentIN = theFastLightVector->
    +                    GetProperty();

dengzy commented 11 years ago

    @@ -780,7 +784,7 @@ void DsG4Scintillation::BuildThePhysicsTable()
                         // Integral pair

                         G4double currentPM = theFastLightVector->
    -                        Energy(0);
    +                        GetPhotonEnergy();

                         G4double currentCII = 0.0;

    @@ -796,11 +800,12 @@ void DsG4Scintillation::BuildThePhysicsTable()
                         // loop over all (photon energy, intensity)
                         // pairs stored for this material

    -                    for (G4int ii=1; ii < theFastLightVector->GetVectorLength(); ++ii) {
    +                    while(++(*theFastLightVector)) {
                             currentPM = theFastLightVector->
    -                            Energy(ii);
    +                            GetPhotonEnergy();

    -                        currentIN=(*theFastLightVector)[ii];
    +                        currentIN=theFastLightVector->  
    +                            GetProperty();

                             currentCII = 0.5 * (prevIN + currentIN);

dengzy commented 11 years ago

低能的迁移代码：

    diff --git a/DetSim/PhysiSim/src/DsPhysConsEM.cc b/DetSim/PhysiSim/src/DsPhysConsEM.cc
    index 70833a7..b629c47 100644
    --- a/DetSim/PhysiSim/src/DsPhysConsEM.cc
    +++ b/DetSim/PhysiSim/src/DsPhysConsEM.cc
    @@ -16,43 +16,16 @@
     #include "G4MuonMinusCaptureAtRest.hh"

     /////////////////////////////// Low Energy EM Process
    -//
    -// New: 2012.12.10
    -// Refer to https://twiki.cern.ch/twiki/bin/view/Geant4/LoweMigration

    -#include "G4PhysicsListHelper.hh"
    +#include "G4LowEnergyRayleigh.hh"
    +#include "G4LowEnergyPhotoElectric.hh"
    +#include "G4LowEnergyCompton.hh"  
    +#include "G4LowEnergyGammaConversion.hh"

    -// gamma
    -#include "G4PhotoElectricEffect.hh"
    -#include "G4LivermorePhotoElectricModel.hh"
    -
    -#include "G4ComptonScattering.hh"
    -#include "G4LivermoreComptonModel.hh"
    -
    -#include "G4GammaConversion.hh"
    -#include "G4LivermoreGammaConversionModel.hh"
    -
    -#include "G4RayleighScattering.hh" 
    -#include "G4LivermoreRayleighModel.hh"
    -// end gamma
    -
    -// e-
    -#include "G4eMultipleScattering.hh"
    -#include "G4UniversalFluctuation.hh"
    -
    -#include "G4eIonisation.hh"
    -#include "G4LivermoreIonisationModel.hh"

dengzy commented 11 years ago

    -// Low Energy Ion
    -// Refer to:
    -//      https://twiki.cern.ch/twiki/bin/view/Geant4/LoweIonParameterized
    -#include "G4ionIonisation.hh"
    -#include "G4IonParametrisedLossModel.hh"
    +#include "G4LowEnergyIonisation.hh"
    +#include "G4LowEnergyBremsstrahlung.hh"

    +#include "G4hLowEnergyIonisation.hh"

dengzy commented 11 years ago

    -//    //fluorescence apply specific cut for fluorescence from photons, electrons
    -//    //and bremsstrahlung photons:
    -//    lowePhot->SetCutForLowEnSecPhotons(m_fluorCut);
    -//    loweIon->SetCutForLowEnSecPhotons(m_fluorCut);
    -//    loweBrem->SetCutForLowEnSecPhotons(m_fluorCut);
    +    G4LowEnergyPhotoElectric* lowePhot = new G4LowEnergyPhotoElectric();
    +    G4LowEnergyIonisation* loweIon  = new G4LowEnergyIonisation();
    +    G4LowEnergyBremsstrahlung* loweBrem = new G4LowEnergyBremsstrahlung();
    +
    +    // note LowEIon uses proton as basis for its data-base, therefore
    +    // cannot specify different LowEnergyIonisation models for different
    +    // particles, but can change model globally for Ion, Alpha and Proton.
    +   
    +    //fluorescence apply specific cut for fluorescence from photons, electrons
    +    //and bremsstrahlung photons:
    +    lowePhot->SetCutForLowEnSecPhotons(m_fluorCut);
    +    loweIon->SetCutForLowEnSecPhotons(m_fluorCut);
    +    loweBrem->SetCutForLowEnSecPhotons(m_fluorCut);

dengzy commented 11 years ago

    @@ -156,17 +104,11 @@ void DsPhysConsEM::ConstructProcess()
                       (particle->GetParticleName() != "chargedgeantino")) {
                 // all other charged particles except geantino
                 G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
    +            G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation();
                 pmanager->AddProcess(aMultipleScattering,-1,1,1);
    -
    -            // Low Energy Ion.
    -            G4ionIonisation* ionIoniProcess = new G4ionIonisation();
    -            G4IonParametrisedLossModel* ionModel = new G4IonParametrisedLossModel();
    -            ionIoniProcess -> SetEmModel(ionModel);
    -            pmanager -> AddProcess(ionIoniProcess, -1, 2, 2) ;
    -
    -            // TODO
    +            pmanager->AddProcess(ahadronLowEIon,-1,2,2);
                 //fluorescence switch off for hadrons :
    -            //ahadronLowEIon->SetFluorescence(false);
    +            ahadronLowEIon->SetFluorescence(false);
             }
         }

dengzy commented 11 years ago

    diff --git a/DetSim/PhysiSim/src/DsPhysConsHadron.cc b/DetSim/PhysiSim/src/DsPhysConsHadron.cc
    index 8031a3c..1f8945b 100644
    --- a/DetSim/PhysiSim/src/DsPhysConsHadron.cc
    +++ b/DetSim/PhysiSim/src/DsPhysConsHadron.cc
    @@ -45,7 +45,7 @@
     #include "G4LFission.hh"
     #include "G4LCapture.hh"

    -#include "G4WHadronElasticProcess.hh"
    +#include "G4UHadronElasticProcess.hh"
     #include "G4HadronElastic.hh"
     #include "G4PreCompoundModel.hh"
     #include "G4BinaryCascade.hh"
    @@ -128,12 +128,6 @@
     #include "G4ProtonInelasticCrossSection.hh"
     #include "G4NeutronInelasticCrossSection.hh"

    -// Hadronic Process
    -//  De-excitation
    -#include "G4Evaporation.hh"
    -#include "G4FermiBreakUp.hh"
    -#include "G4StatMF.hh"
    -
     DsPhysConsHadron::DsPhysConsHadron(const G4String& name): G4VPhysicsConstructor(name)
     {
         m_minEnergyForMultiFrag=3.0*MeV;
    @@ -473,7 +467,7 @@ void DsPhysConsHadron::ConstructProcess()
         pmanager = G4Neutron::Neutron()->GetProcessManager();

         // elastic scattering
    -    G4WHadronElasticProcess* theHadronElasticProcess = new G4WHadronElasticProcess("neutronElastic");
    +    G4UHadronElasticProcess* theHadronElasticProcess = new G4UHadronElasticProcess("neutronElastic");
         G4HadronElastic* theNeutronElasticModel = new G4HadronElastic();
         G4NeutronHPElastic* theNeutronHPElastic = new G4NeutronHPElastic();
         //theNeutronHPElastic->SetMinEnergy(0.*MeV);

mirguest commented 11 years ago

很有可能与强子过程有关。

mirguest commented 11 years ago

需要研究，中子的Elastic过程应该选用哪个比较好。

dengzy commented 11 years ago

其实，是质子的Low Energy Elastic进入了死循环。原因在于，对低能的代码进行迁移时，考虑了低能的电离过程。而现在，把质子的电离过程单独提取出来，即可。