Ignore:
Timestamp:
Nov 25, 2009, 5:13:58 PM (16 years ago)
Author:
garnier
Message:

update CVS release candidate geant4.9.3.01

Location:
trunk/source/processes/electromagnetic/adjoint
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/adjoint/GNUmakefile

    r966 r1196  
    1 # $Id: GNUmakefile,v 1.2 2008/11/14 20:47:47 vnivanch Exp $
     1# $Id: GNUmakefile,v 1.4 2009/11/11 11:23:45 vnivanch Exp $
    22# --------------------------------------------------------------------
    33# GNUmakefile for electromagnetic sub-library.  G.Cosmo, 14/11/2008.
    … …  
    1818            -I$(G4BASE)/geometry/management/include \
    1919            -I$(G4BASE)/geometry/volumes/include \
     20            -I$(G4BASE)/geometry/navigation/include \
    2021            -I$(G4BASE)/track/include \
    2122            -I$(G4BASE)/processes/management/include \

  • trunk/source/processes/electromagnetic/adjoint/History

    r966 r1196  
    1 $Id: History,v 1.1 2008/11/14 19:54:40 gcosmo Exp $
     1$Id: History,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1818     ----------------------------------------------------------
    1919
     2020 Nov 2009: L.Desorgher (emadjoint-V09-02-01)
     21-Replace C++ type by G4 type where needed and adding of G4disclaimer
     22
     23
     2410 Nov 2009: L.Desorgher (emadjoint-V09-02-00)
     25- Commit of the electromagnetic  adjoint processes for the release of the all adjoint machinery into Geant4.
     26  Compared to the first commit, all e- processes have been improved and adjoint proton and ion ionisation have been added.
     27  The use of adjoint cross section matrices can be now limited only to e- Ionisation and Ion ionisation.
     28  The GNUmakefile has been modified by adding  -I$(G4BASE)/geometry/navigation/include  in CPPFLAGS.       
     29
    203014 Nov 2008: G.Cosmo (emadjoint-V09-01-00)
    2131- First commit.

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointAlongStepWeightCorrection.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointAlongStepWeightCorrection.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointAlongStepWeightCorrection.hh
     30//      Class:          G4AdjointAlongStepWeightCorrection
    2831//      Author:         L. Desorgher
    29 //      Date:           10 May 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              10 May 2007 creation by L. Desorgher           
     40//              10 May 2007 creation by L. Desorgher 
     41//              October 2009 implementation of the mode where the total adjoint and forward cross sections are equivalent. L. Desorgher         
    3842//
    3943//-------------------------------------------------------------
    4044//      Documentation:
    4145//              Continuous processes acting on adjoint particles to correct continuously their weight during the adjoint reverse tracking.
     46//              Thi process is needed whene the adjoint cross section are not scaled such that the total adjoint cross section match the total forward cross section.
     47//              By default the mode where the total adjoint cross section is equal to the total forward cross section is used an therefore this along step weight
     48//              correction factor is 1.
     49//              However in some cases (some energy ranges) the total forward cross section or the total adjoint cross section can be null, in this case the along step
     50//              weight correction is neede and is given by exp(-(Sigma_tot_adj-Sigma_tot_fwd).dx)
     51//             
     52//
    4253//
    4354
    … …  
    120131    currentMaterial = couple->GetMaterial();
    121132    currentMaterialIndex = couple->GetIndex();
    122     //G4cout<<"Define Material"<<std::endl;
    123     //if(!meanFreePath) ResetNumberOfInteractionLengthLeft();
     133 
    124134  }
    125135}
    126136
    127137
    128 ///////////////////////////////////////////////////////
    129 //
    130 inline G4double G4AdjointAlongStepWeightCorrection::GetContinuousStepLimit(const G4Track& track,
    131                 G4double , G4double , G4double& )
    132 {
    133   G4double x = DBL_MAX;
    134   DefineMaterial(track.GetMaterialCutsCouple());
    135   preStepKinEnergy = track.GetKineticEnergy();
    136   return x;
    137 }
     138
    138139#endif

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointBremsstrahlungModel.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointBremsstrahlungModel.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointBremsstrahlungModel.hh
     30//      Class:          G4AdjointBremsstrahlungModel
    2831//      Author:         L. Desorgher
    29 //      Date:           15 June 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              15 June 2007 creation by L. Desorgher. Adapted from G4eBremsstrahlungModel             
     40//              15 June 2007 creation by L. Desorgher. Adapted from G4eBremsstrahlungModel 
     41//              20-10-2009 Remove all the screening effect that are not considered in the direct models blow 10 GeV. L.Desorgher
     42//              4-11-2009  Implement the use of a simple biased differential cross section (C(Z)/Egamma) allowing a rapid computation of adjoint CS
     43//                         and rapid sampling of adjoint secondaries. By this way cross section matrices are not used anymore, avoiding a rather
     44//                         time consuming computation of adjoint brem cross section matrices for each material at initialisation. This mode is switch on/off
     45//                         by selecting SetUseMatrix(false)/ SetUseMatrix(true) in the constructor. L.Desorgher
     46//                                     
    3847//
    3948//-------------------------------------------------------------
    … …  
    4958#include "G4VEmAdjointModel.hh"
    5059#include "G4eBremsstrahlungModel.hh"
     60//#include "G4PenelopeBremsstrahlungModel.hh"
     61#include "G4PhysicsTable.hh"
     62//#include "G4EmModelManager.hh"
    5163class G4Timer;
    5264class G4AdjointBremsstrahlungModel: public G4VEmAdjointModel
    … …  
    6072                                G4bool IsScatProjToProjCase,
    6173                                G4ParticleChange* fParticleChange);
    62  
     74  void RapidSampleSecondaries(const G4Track& aTrack,
     75                                G4bool IsScatProjToProjCase,
     76                                G4ParticleChange* fParticleChange);
    6377  virtual G4double DiffCrossSectionPerVolumePrimToSecond(
    6478                                      const G4Material* aMaterial,
    … …  
    6680                                      G4double kinEnergyProd // kinetic energy of the secondary particle
    6781                                      );
    68   G4double DiffCrossSectionPerVolumePrimToSecond1(
     82  G4double DiffCrossSectionPerVolumePrimToSecondApproximated1(
    6983                                      const G4Material* aMaterial,
    7084                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    7185                                      G4double kinEnergyProd // kinetic energy of the secondary particle
    7286                                      );
    73   G4double DiffCrossSectionPerVolumePrimToSecond2(
     87  G4double DiffCrossSectionPerVolumePrimToSecondApproximated2(
    7488                                      const G4Material* aMaterial,
    7589                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    7690                                      G4double kinEnergyProd // kinetic energy of the secondary particle
    7791                                      );
    78   G4double DiffCrossSectionPerVolumePrimToSecond3(
    79                                       const G4Material* aMaterial,
    80                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    81                                       G4double kinEnergyProd // kinetic energy of the secondary particle
    82                                       );                                     
    83   void DefineDirectBremModel(G4eBremsstrahlungModel* aModel);
    84   inline void SetdCSModel(G4String aString) {ModeldCS=aString;}
    85                                      
     92  virtual G4double AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
     93                                             G4double primEnergy,
     94                                             G4bool IsScatProjToProjCase);
    8695 
    87 private:
    88  
    89   void InitialiseParameters();
    90   G4double SupressionFunction(const G4Material* material, G4double tkin,
    91                                     G4double gammaEnergy);
    92                                    
     96 
     97 // private void InitialiseFwdModels();
     98
    9399
    94100private: 
    95   G4eBremsstrahlungModel* theDirectBremModel;
     101   G4eBremsstrahlungModel* theDirectStdBremModel;
     102  //G4PenelopeBremsstrahlungModel* theDirectPenelopeBremModel;
     103
    96104
    97105  G4double highKinEnergy;
    98106  G4double lowKinEnergy;
    99   G4double probsup;
    100   G4double MigdalConstant;
    101   G4double LPMconstant;
    102   G4double highEnergyTh;
    103   G4bool   theLPMflag;
    104   G4bool   isElectron;
     107 
    105108 
    106   //Vector
     109  std::vector<G4DataVector*> partialSumSigma;
    107110 
    108   std::vector<float> FZ;
    109   std::vector<float> ah1;
    110   std::vector<float> ah2;
    111   std::vector<float> ah3;
    112  
    113   std::vector<float> bh1;
    114   std::vector<float> bh2;
    115   std::vector<float> bh3;
    116  
    117   std::vector<float> al0;
    118   std::vector<float> al1;
    119   std::vector<float> al2;
    120  
    121   std::vector<float> bl0;
    122   std::vector<float> bl1;
    123   std::vector<float> bl2;
     111 
    124112 
    125113  std::vector<float> SigmaPerAtom;
    126114  G4Timer* theTimer;
    127115 
    128   G4String ModeldCS;
     116  G4double MigdalConstant;
     117  G4double lastCZ;
     118 
     119
     120 /*
     121  G4bool UsePenelopeModel;
     122  G4EmModelManager* theEmModelManagerForFwdModels;
     123  G4bool isPenelopeModelInitialised ;
     124 */
    129125 
    130126 
    131127};
    132128
     129
    133130#endif

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointCSManager.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointCSManager.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointCSManager.hh
     30//      Class:          G4AdjointCSManager
    2831//      Author:         L. Desorgher
    29 //      Date:           1st April 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3639//      ChangeHistory:
    3740//              1st April 2007 creation by L. Desorgher                 
     41//             
     42//              September-October 2009. Implementation of the mode where the adjoint cross sections are scaled such that the total used adjoint cross sections is in
     43//              most of the cases equal to the total forward cross section. L.Desorgher
    3844//
    3945//-------------------------------------------------------------
    4046//      Documentation:
    4147//              Is responsible for the management of all adjoint cross sections matrices, and for the computation of the total forward and adjoint cross sections.
    42 //              Total adjoint and forward cross sections are needed to correct continuously the weight of a particle after a tracking step.
     48//              Total adjoint and forward cross sections are needed to correct the weight of a particle after a tracking step or after the occurence of a reverse reaction.
    4349//              It is also used to sample an adjoint secondary from a given adjoint cross section matrix.
    4450//
    … …  
    8288        void RegisterAdjointParticle(G4ParticleDefinition* aPartDef);
    8389       
    84         //Building of thr CS Matrices and Total Forward and Adjoint LambdaTables
     90        //Building of the CS Matrices and Total Forward and Adjoint LambdaTables
    8591        //----------------------------------------------------------------------
    8692       
    … …  
    8995       
    9096       
    91         //Get TotalCrossSections form Total Lambda Tables
     97        //Get TotalCrossSections form Total Lambda Tables, Needed for Weight correction and scaling of the
    9298        //-------------------------------------------------
    9399        G4double GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
    94100                                     const G4MaterialCutsCouple* aCouple);
    95101        G4double GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
    96                                      const G4MaterialCutsCouple* aCouple);                           
    97        
     102                                     const G4MaterialCutsCouple* aCouple);
     103       
     104       
     105        void GetEminForTotalCS(G4ParticleDefinition* aPartDef,
     106                                     const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd);
     107        void GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
     108                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max);
     109        void GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
     110                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max);
     111                                     
     112                                                                                     
     113       
     114        //CrossSection Correction  1 or FwdCS/AdjCS following the G4boolean value of forward_CS_is_used and forward_CS_mode
     115        //-------------------------------------------------
     116        G4double GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used, G4double& fwd_TotCS);
     117       
     118       
     119        //Cross section mode
     120        //------------------
     121        inline void SetFwdCrossSectionMode(G4bool aBool){forward_CS_mode=aBool;}
    98122       
    99123       
    … …  
    103127        G4double GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
    104128                                     const G4MaterialCutsCouple* aCouple, G4double step_length);
    105         G4double GetPostStepWeightCorrection(G4ParticleDefinition* aPrimPartDef, G4ParticleDefinition* aSecondPartDef,
    106                                             G4double EkinPrim,G4double EkinSecond,
    107                                      const G4MaterialCutsCouple* aCouple);
     129        G4double GetPostStepWeightCorrection();
    108130                                     
    109131       
    110         double ComputeAdjointCS(G4Material* aMaterial,
     132       
     133       
     134        //Method Called by the adjoint model to get there CS, if not precised otherwise
     135        //-------------------------------
     136       
     137        G4double ComputeAdjointCS(G4Material* aMaterial,
    111138                                            G4VEmAdjointModel* aModel,
    112139                                            G4double PrimEnergy,
    113140                                            G4double Tcut,
    114141                                            G4bool IsScatProjToProjCase,
    115                                             std::vector<double>&
     142                                            std::vector<G4double>&
    116143                                                 AdjointCS_for_each_element);
    117144                                         
    118        
     145        //Method Called by the adjoint model to sample the secondary energy form the CS matrix
     146        //--------------------------------------------------------------------------------
    119147        G4Element*  SampleElementFromCSMatrices(G4Material* aMaterial,
    120148                                                       G4VEmAdjointModel* aModel,
    … …  
    122150                                                        G4double Tcut,
    123151                                                       G4bool IsScatProjToProjCase);
    124         G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple* aMatCutCouple,G4ParticleDefinition* aPart,G4double PrimEnergy);                                                     
     152                                                       
     153                                                       
     154        //Total Adjoint CS  is computed at initialisation phase
     155        //-----------------------------------------------------                                       
     156        G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple* aMatCutCouple,G4ParticleDefinition* aPart,G4double PrimEnergy);
     157       
     158       
     159       
     160                                                               
    125161        G4ParticleDefinition* GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef);
    126162        G4ParticleDefinition* GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef);
    … …  
    130166        inline void SetTmax(G4double aVal){Tmax=aVal;}
    131167        inline void SetNbins(G4int aInt){nbins=aInt;}
    132        
    133        
    134        
    135         //inline
    136         inline void ConsiderContinuousWeightCorrection(G4bool aBool){consider_continuous_weight_correction=aBool;}
    137         inline void ConsiderPoststepWeightCorrection(G4bool aBool){consider_poststep_weight_correction=aBool;}
    138        
    139        
     168        inline void SetIon(G4ParticleDefinition* adjIon,
     169                           G4ParticleDefinition* fwdIon) {theAdjIon=adjIon; theFwdIon =fwdIon;}
    140170       
    141171       
    … …  
    163193        std::vector< size_t> listOfIndexOfAdjointEMModelInAction;
    164194        std::vector< G4bool> listOfIsScatProjToProjCase;
    165         std::vector< std::vector<double> > lastAdjointCSVsModelsAndElements;
     195        std::vector< std::vector<G4double> > lastAdjointCSVsModelsAndElements;
    166196        G4bool CrossSectionMatrixesAreBuilt;
     197        size_t  currentParticleIndex;
     198        G4ParticleDefinition* currentParticleDef;
    167199       
    168200        //total adjoint and total forward cross section table in function of material and in function of adjoint particle type
    … …  
    170202        std::vector<G4PhysicsTable*>        theTotalForwardSigmaTableVector;
    171203        std::vector<G4PhysicsTable*>        theTotalAdjointSigmaTableVector;
     204        std::vector< std::vector<G4double> >    EminForFwdSigmaTables; 
     205        std::vector< std::vector<G4double> >    EminForAdjSigmaTables; 
     206        std::vector< std::vector<G4double> >    EkinofFwdSigmaMax;     
     207        std::vector< std::vector<G4double> >    EkinofAdjSigmaMax;
     208       
     209       
     210       
    172211         
    173212        //list of forward G4VEMLossProcess and of G4VEMProcess for the different adjoint particle
    … …  
    177216       
    178217        //list of adjoint particles considered
    179        
     218        //--------------------------------------------------------------
    180219        std::vector< G4ParticleDefinition*> theListOfAdjointParticlesInAction;
    181220       
    … …  
    191230        size_t  currentMatIndex;
    192231       
    193         int verbose;
    194        
    195        
    196         //Weight correction
    197         //------------------
    198         G4bool consider_continuous_weight_correction;
    199         G4bool consider_poststep_weight_correction;
     232        G4int verbose;
     233       
     234       
     235       
     236       
     237        //Two CS mode are possible :forward_CS_mode = false the Adjoint CS are used as it is implying a AlongStep Weight Correction.
     238        //                         :forward_CS_mode = true the Adjoint CS are scaled to have the total adjoint CS eual to the fwd one implying a PostStep Weight Correction.
     239        //                                           For energy range  where the total FwdCS or the total adjoint CS are null, the scaling is not possble and
     240        //                                           forward_CS_is_used is set to false         
     241        //--------------------------------------------
     242        G4bool forward_CS_is_used;
     243        G4bool forward_CS_mode;
     244       
     245        //Adj and Fwd CS values for re-use
     246        //------------------------
     247       
     248        G4double PreadjCS,PostadjCS;
     249        G4double PrefwdCS,PostfwdCS;
     250        G4double LastEkinForCS;
     251        G4double LastCSCorrectionFactor;
     252        G4ParticleDefinition* lastPartDefForCS;
     253       
     254        //Ion
     255        //----------------
     256        G4ParticleDefinition* theAdjIon; //at the moment Only one ion can be considered by simulation
     257        G4ParticleDefinition* theFwdIon;
     258        G4double massRatio;
     259       
     260       
     261       
    200262
    201263  private:
    202264        G4AdjointCSManager(); 
    203265        void DefineCurrentMaterial(const G4MaterialCutsCouple* couple);
    204         double ComputeAdjointCS(G4double aPrimEnergy, G4AdjointCSMatrix* anAdjointCSMatrix, G4double Tcut);
     266        void DefineCurrentParticle(const G4ParticleDefinition* aPartDef);
     267        G4double ComputeAdjointCS(G4double aPrimEnergy, G4AdjointCSMatrix* anAdjointCSMatrix, G4double Tcut);
    205268
    206269};

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointCSMatrix.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointCSMatrix.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointCSMatrix.hh
     30//      Class:          G4AdjointCSMatrix.hh
    2831//      Author:         L. Desorgher
    29 //      Date:           1st April 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    6568        ~G4AdjointCSMatrix();
    6669
    67         ////////////
    68         // Methods
    69         ////////////
     70        //////////////
     71        // Methods  //
     72        //////////////
    7073        void Clear();
    71         void AddData(G4double aPrimEnergy,G4double aCS, std::vector< G4double>* aLogSecondEnergyVector,
    72                                                         std::vector< G4double>* aLogProbVector,size_t n_pro_decade=0); 
     74        void AddData(G4double aPrimEnergy,G4double aCS, std::vector< double>* aLogSecondEnergyVector,
     75                                                        std::vector< double>* aLogProbVector,size_t n_pro_decade=0);   
    7376       
    74         bool GetData(unsigned int i, G4double& aPrimEnergy,G4double& aCS,G4double& log0, std::vector< G4double>*& aLogSecondEnergyVector,
    75                                                                       std::vector< G4double>*& aLogProbVector,
     77        G4bool GetData(unsigned int i, G4double& aPrimEnergy,G4double& aCS,G4double& log0, std::vector< double>*& aLogSecondEnergyVector,
     78                                                                      std::vector< double>*& aLogProbVector,
    7679                                                                      std::vector< size_t>*& aLogProbVectorIndex);
    7780       
    78         inline std::vector< G4double >* GetLogPrimEnergyVector(){return &theLogPrimEnergyVector;}
    79         inline std::vector< G4double >* GetLogCrossSectionvector(){return &theLogCrossSectionVector;}
     81        inline std::vector< double>* GetLogPrimEnergyVector(){return &theLogPrimEnergyVector;}
     82        inline std::vector< double>* GetLogCrossSectionvector(){return &theLogCrossSectionVector;}
    8083        inline G4double GetDlog(){return dlog;}         
    8184        inline G4bool IsScatProjToProjCase(){return is_scat_proj_to_proj_case;}
    … …  
    8790        // we did first try to use G4PhysicsOrderedVector but they are not general enough for our purpose
    8891       
    89         std::vector< G4double > theLogPrimEnergyVector;
    90         std::vector< G4double > theLogCrossSectionVector; //Adjoint Cross sections in function of primary energy
    91         std::vector< std::vector< G4double >* > theLogSecondEnergyMatrix;
    92         std::vector< std::vector< G4double >* > theLogProbMatrix; //Each column represents the integrated probability of getting a secondary
     92        std::vector< double> theLogPrimEnergyVector;
     93        std::vector< double> theLogCrossSectionVector; //Adjoint Cross sections in function of primary energy
     94        std::vector< std::vector< double>* > theLogSecondEnergyMatrix;
     95        std::vector< std::vector< double>* > theLogProbMatrix; //Each column represents the integrated probability of getting a secondary
    9396                                                                      // in function of their energy
    94         std::vector< std::vector< size_t >* > theLogProbMatrixIndex; //index of euqidistant LogProb
    95         std::vector< G4double > log0Vector;
     97        std::vector< std::vector< size_t >* > theLogProbMatrixIndex; //index of equidistant LogProb
     98        std::vector< double> log0Vector;
    9699       
    97100        unsigned int nb_of_PrimEnergy;

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointComptonModel.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointComptonModel.hh,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointComptonModel.hh
     30//      Class:          G4AdjointComptonModel
    2831//      Author:         L. Desorgher
    29 //      Date:           1 September 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              1 September 2007 creation by L. Desorgher               
     40//              1 September 2007 creation by L. Desorgher
     41//              11 November 2009 Implement the use of approximated diffCS as an alternative of CSMatrix.               
    3842//
    3943//-------------------------------------------------------------
    4044//      Documentation:
    41 //              Model for the adjoint compton scattering
     45//              Model for the adjoint compton scattering.
    4246//
    4347
    … …  
    4852#include "globals.hh"
    4953#include "G4VEmAdjointModel.hh"
     54#include "G4VEmProcess.hh"
    5055class G4AdjointComptonModel: public G4VEmAdjointModel
    5156
    … …  
    5863 
    5964  virtual void SampleSecondaries(const G4Track& aTrack,
     65                                G4bool IsScatProjToProjCase,
     66                                G4ParticleChange* fParticleChange);
     67  void RapidSampleSecondaries(const G4Track& aTrack,
    6068                                G4bool IsScatProjToProjCase,
    6169                                G4ParticleChange* fParticleChange);
    … …  
    7179                                      G4double Z,
    7280                                      G4double A = 0.);
     81 
    7382  virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy);
    7483  virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy);
    75 
     84 
     85 
     86  virtual G4double AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
     87                                             G4double primEnergy,
     88                                             G4bool IsScatProjToProjCase);
    7689
    7790 
    7891 
    79                              
     92  inline void SetDirectProcess(G4VEmProcess* aProcess){theDirectEMProcess = aProcess;};                       
    8093 
    8194private:
     95  G4VEmProcess* theDirectEMProcess;
     96  G4double G4direct_CS;
    8297 
    8398   

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointInterpolator.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointInterpolator.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointInterpolator.hh
     30//      Module:         G4AdjointInterpolator
    2831//      Author:         L. Desorgher
    29 //      Date:           1st April 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    7780       
    7881       
    79         size_t FindPosition(G4double& x,std::vector<double>& x_vec,size_t ind_min=0, size_t ind_max=0);
     82        size_t FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t ind_min=0, size_t ind_max=0);
    8083       
    81         size_t FindPositionForLogVector(G4double& x,std::vector<double>& x_vec);
     84        size_t FindPositionForLogVector(G4double& x,std::vector<G4double>& x_vec);
    8285       
    83         G4double Interpolate(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,G4String InterPolMethod="Log"); //xvec should monotically increase                       
     86        G4double Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod="Log"); //xvec should monotically increase                   
    8487       
    85         G4double InterpolateWithIndexVector(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,
     88        G4double InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
    8689                                            std::vector<size_t>& index_vec, G4double x0,G4double dx); //xvec should monotically increase                       
    8790       
    8891       
    89         G4double InterpolateForLogVector(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec);
     92        G4double InterpolateForLogVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec);
    9093   
    9194   private:

  • trunk/source/processes/electromagnetic/adjoint/include/G4AdjointPhotoElectricModel.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointPhotoElectricModel.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointPhotoElectricModel.hh
     30//      Module:         G4AdjointPhotoElectricModel
    2831//      Author:         L. Desorgher
    29 //      Date:           10 October 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              1 September 2007 creation by L. Desorgher               
     40//              -1 September 2007 creation by L. Desorgher 
     41//             
     42//              -January 2009. L. Desorgher     
     43//               Put a higher limit on the CS to avoid a high rate of  Inverse Photo e- effect at low energy. The very high adjoint CS of the reverse
     44//               photo electric reaction produce a high rate of reverse photo electric reaction in the inner side of a shielding for eaxmple, the correction of this occurence
     45//               by weight correction in the StepDoIt method is not statistically sufficient at small energy. The problem is partially solved by setting an higher CS limit
     46//               and compensating it by an extra weight correction factor. However when  coupling it with other reverse  processes the reverse photo-electric is still
     47//               the source of very occasional high weight that decrease the efficiency of the computation. A way to solve this problemn is still needed but is difficult
     48//               to find as it happens in rarea case but does give a weighrt that is outside the noemal distribution. (Very Tricky!) 
     49//             
     50//              -October 2009   Correction of Element sampling. L. Desorgher           
    3851//
    3952//-------------------------------------------------------------
    … …  
    6174                                G4bool IsScatProjToProjCase,
    6275                                G4ParticleChange* fParticleChange);
    63  
    6476  virtual G4double AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
    6577                                G4double primEnergy,
    … …  
    7385                                                       DefineDirectEMModel(aModel);}                                 
    7486 
     87  virtual void CorrectPostStepWeight(G4ParticleChange* fParticleChange,
     88                                     G4double old_weight,
     89                                     G4double adjointPrimKinEnergy,
     90                                     G4double projectileKinEnergy,
     91                                     G4bool IsScatProjToProjCase);
    7592 
    7693 
    … …  
    7895  G4double  xsec[40];
    7996  G4double  totAdjointCS;
     97  G4double  totBiasedAdjointCS;
     98  G4double  factorCSBiasing;
     99  G4double  pre_step_AdjointCS;
     100  G4double  post_step_AdjointCS;
     101 
     102 
    80103  G4double  shell_prob[40][40];
    81104 

  • trunk/source/processes/electromagnetic/adjoint/include/G4ContinuousGainOfEnergy.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4ContinuousGainOfEnergy.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4ContinuousGainOfEnergy.hh
     30//      Class:          G4ContinuousGainOfEnergy
    2831//      Author:         L. Desorgher
    29 //      Date:           10 May 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              10 May 2007 creation by L. Desorgher           
     40//              -10 May 2007 creation by L. Desorgher 
     41//              -February-March 2009 Update for protons by L.Desorgher
     42//              -July August 2009  Update for ion by L.Desorgher
    3843//
    3944//-------------------------------------------------------------
    4045//      Documentation:
    41 //              Continuous process acting on adjoint particles to compute the continuous gain of energy of charged particels whern they are tracked back!
     46//              Continuous process acting on adjoint particles to compute the continuous gain of energy of charged particles when they are tracked back!
     47//             
    4248//
    4349#ifndef G4ContinuousGainOfEnergy_h
    … …  
    5258#include "G4ParticleChange.hh"
    5359#include "G4VEnergyLossProcess.hh"
     60#include "G4ProductionCutsTable.hh"
    5461
    5562
    … …  
    106113  inline void SetDirectEnergyLossProcess(G4VEnergyLossProcess* aProcess){theDirectEnergyLossProcess=aProcess;}; 
    107114 
    108   inline void SetDirectParticle(G4ParticleDefinition* p){theDirectPartDef=p;};
     115  void SetDirectParticle(G4ParticleDefinition* p);
    109116
    110117protected:
    … …  
    116123
    117124  void DefineMaterial(const G4MaterialCutsCouple* couple);
     125  void SetDynamicMassCharge(const G4Track& track, G4double energy);
     126 
    118127 
    119128  // hide  assignment operator
    … …  
    128137  const G4MaterialCutsCouple* currentCouple;
    129138  size_t   currentMaterialIndex;
     139  size_t   currentCoupleIndex;
    130140  G4double currentTcut;
     141  G4double currentCutInRange;
    131142  G4double preStepKinEnergy;
     143 
    132144 
    133145 
    … …  
    142154  G4bool is_integral;
    143155 
     156  //adding for Ions
     157  //----------------
     158  G4bool IsIon;
     159  G4double massRatio;
     160  G4double chargeSqRatio;
     161  G4VEmModel* currentModel;
     162  G4double preStepChargeSqRatio;
     163  G4double preStepScaledKinEnergy;
     164
     165 
     166 
     167 
    144168 
    145169};
    … …  
    153177    currentCouple   = couple;
    154178    currentMaterial = couple->GetMaterial();
    155     currentMaterialIndex = couple->GetIndex();
    156     currentTcut = couple->GetProductionCuts()->GetProductionCut(theDirectPartDef->GetParticleName());
    157     //G4cout<<"Define Material"<<std::endl;
     179    currentCoupleIndex = couple->GetIndex();
     180    currentMaterialIndex = currentMaterial->GetIndex();
     181   
     182    size_t idx=1;
     183    const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
     184    currentTcut=(*aVec)[currentCoupleIndex];
     185    currentCutInRange = couple->GetProductionCuts()->GetProductionCut(theDirectPartDef->GetParticleName());
     186    //G4cout<<"Define Material"<<G4endl;
    158187    //if(!meanFreePath) ResetNumberOfInteractionLengthLeft();
    159188  }
    160189}
    161 ///////////////////////////////////////////////////////
    162 //
    163 inline G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
    164                 G4double , G4double , G4double& )
    165 {
    166   G4double x = DBL_MAX;
    167   x=.1*mm;
    168  
    169   //G4cout<<x<<std::endl;
    170   DefineMaterial(track.GetMaterialCutsCouple());
    171   preStepKinEnergy = track.GetKineticEnergy();
    172   G4double maxE=1.2*preStepKinEnergy;
    173   G4double r = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
    174   G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
    175   x=std::max(r1-r,.1);
    176  
    177   return x;
    178  
    179  
    180 }
     190
    181191#endif

  • trunk/source/processes/electromagnetic/adjoint/include/G4InversePEEffect.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4InversePEEffect.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4AdjointPEEffect.hh
     30//      Module:         G4InversePEEffect.hh
    2831//      Author:         L. Desorgher
    29 //      Date:           25 October 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    4548#define G4InversePEEffect_h 1
    4649
    47 #include "G4VAdjointInverseScattering.hh"
     50#include "G4VAdjointReverseReaction.hh"
    4851#include "globals.hh"
    4952class G4AdjointPhotoElectricModel;
    50 class G4InversePEEffect: public G4VAdjointInverseScattering
     53class G4InversePEEffect: public G4VAdjointReverseReaction
    5154
    5255{

  • trunk/source/processes/electromagnetic/adjoint/include/G4VEmAdjointModel.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4VEmAdjointModel.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    27 //      Module:         G4VEMAdjointModel.hh
     30//      Module:         G4VEMAdjointModel
    2831//      Author:         L. Desorgher
    29 //      Date:           1st April 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3538// --------------
    3639//      ChangeHistory:
    37 //              1st April 2007 creation by L. Desorgher                 
     40//              10 September 2009 Move to a virtual class. L. Desorgher
     41//              1st April 2007 creation by L. Desorgher                 
    3842//
    3943//-------------------------------------------------------------
    4044//      Documentation:
    41 //              Base class for Adjoint model
    42 //
     45//              Base class for Adjoint EM model. It is based on the use of direct G4VEmModel.
     46//
     47
    4348
    4449#ifndef G4VEmAdjointModel_h
    … …  
    6974{
    7075
    71 public:
     76public: // public methods
    7277
    7378  G4VEmAdjointModel(const G4String& nam);
    … …  
    7681
    7782  //------------------------------------------------------------------------
    78   // Virtual methods to be implemented for the concrete model
     83  // Virtual methods to be implemented for the sample secondaries concrete model
    7984  //------------------------------------------------------------------------
    80 
    81   //virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&) = 0;
    82 
    83 
     85 
     86  //virtual void Initialise()=0;
     87 
    8488  virtual void SampleSecondaries(const G4Track& aTrack,
    8589                                G4bool IsScatProjToProjCase,
    86                                 G4ParticleChange* fParticleChange);
    87  
     90                                G4ParticleChange* fParticleChange)=0;
    8891 
    8992
    … …  
    9194  // Methods for adjoint processes; may be overwritten if needed; 
    9295  //------------------------------------------------------------------------
     96 
     97
    9398  virtual G4double AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
    9499                                             G4double primEnergy,
    … …  
    121126                                      );
    122127 
    123                                      
    124  
    125  
    126   G4double DiffCrossSectionFunction1(G4double kinEnergyProj);
    127   G4double DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd);
    128  
    129   G4double DiffCrossSectionFunction2(G4double kinEnergyProj);
    130  
    131   std::vector< std::vector< G4double >* >  ComputeAdjointCrossSectionVectorPerAtomForSecond(     
     128 
     129  //Energy limits of adjoint secondary
     130  //------------------
     131 
     132  virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy);
     133  virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut=0);
     134  virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy);
     135  virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy);
     136 
     137 
     138 
     139  //Other Methods
     140  //---------------
     141 
     142  void  DefineCurrentMaterial(const G4MaterialCutsCouple* couple);
     143 
     144 
     145  std::vector< std::vector< double>* >  ComputeAdjointCrossSectionVectorPerAtomForSecond(     
    132146                                G4double kinEnergyProd,
    133147                                G4double Z,
    … …  
    135149                                G4int nbin_pro_decade=10
    136150                                );
    137   std::vector< std::vector< G4double >* >  ComputeAdjointCrossSectionVectorPerAtomForScatProj(     
     151  std::vector< std::vector< double>* >  ComputeAdjointCrossSectionVectorPerAtomForScatProj(     
    138152                                G4double kinEnergyProd,
    139153                                G4double Z,
    … …  
    142156                                );
    143157 
    144   std::vector< std::vector< G4double >* >  ComputeAdjointCrossSectionVectorPerVolumeForSecond(     
     158  std::vector< std::vector< double>* >  ComputeAdjointCrossSectionVectorPerVolumeForSecond(     
    145159                                G4Material* aMaterial,
    146160                                G4double kinEnergyProd,
    147161                                G4int nbin_pro_decade=10
    148162                                );
    149   std::vector< std::vector< G4double >* >  ComputeAdjointCrossSectionVectorPerVolumeForScatProj(     
     163  std::vector< std::vector< double>* >  ComputeAdjointCrossSectionVectorPerVolumeForScatProj(     
    150164                                G4Material* aMaterial,
    151165                                G4double kinEnergyProd,
    152166                                G4int nbin_pro_decade=10
    153167                                );
    154                                
    155  
    156                                                        
    157   virtual G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double prim_energy,G4bool IsScatProjToProjCase);     
    158   virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase);     
    159   void CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy);     
    160  
    161   //Set/Get methods
    162   //------------------
    163  
    164   virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy);
    165   virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut=0);
    166   virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy);
    167   virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy);
    168  
    169   virtual void SetCSBiasingFactor(G4double aVal) {CS_biasing_factor = aVal;}
    170  
    171  
    172  
    173  
    174  
    175  
    176 public: 
     168 
     169
    177170 
    178171  inline void SetCSMatrices(std::vector< G4AdjointCSMatrix* >* Vec1CSMatrix, std::vector< G4AdjointCSMatrix* >* Vec2CSMatrix){
    … …  
    191184  inline G4double GetLowEnergyLimit(){return LowEnergyLimit;}
    192185 
    193   inline void SetHighEnergyLimit(G4double aVal){HighEnergyLimit=aVal;}
    194  
    195   inline void SetLowEnergyLimit(G4double aVal){LowEnergyLimit=aVal;}
    196  
    197   inline void SetCorrectWeightMode(G4bool aBool){CorrectWeightMode=aBool;};
    198  
    199   inline void SetApplyBiasing(G4bool aBool){ApplyBiasing=aBool;};
    200  
    201                              
     186  void SetHighEnergyLimit(G4double aVal);
     187 
     188  void SetLowEnergyLimit(G4double aVal);
     189 
    202190  inline void DefineDirectEMModel(G4VEmModel* aModel){theDirectEMModel = aModel;}
    203191 
    204   inline void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart){
    205         theAdjEquivOfDirectPrimPartDef=aPart;
    206         if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
    207                                         theDirectPrimaryPartDef=G4Electron::Electron();
    208         if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
    209                                         theDirectPrimaryPartDef=G4Gamma::Gamma();
    210  
    211   }
     192  void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart);
    212193 
    213194  inline void SetAdjointEquivalentOfDirectSecondaryParticleDefinition(G4ParticleDefinition* aPart){
    … …  
    217198  inline void SetSecondPartOfSameType(G4bool aBool){second_part_of_same_type =aBool;}
    218199 
    219   bool GetSecondPartOfSameType(){return second_part_of_same_type;}
     200  inline G4bool GetSecondPartOfSameType(){return second_part_of_same_type;}
    220201 
    221202  inline void SetUseMatrix(G4bool aBool) { UseMatrix = aBool;}
    … …  
    225206 
    226207  inline void SetApplyCutInRange(G4bool aBool){ ApplyCutInRange = aBool;}
    227   inline void SetIsIonisation(G4bool aBool){ IsIonisation = aBool;}
    228  
    229208  inline G4bool GetUseMatrix() {return UseMatrix;}
    230209  inline G4bool GetUseMatrixPerElement(){ return UseMatrixPerElement;}
    231210  inline G4bool GetUseOnlyOneMatrixForAllElements(){ return UseOnlyOneMatrixForAllElements;}
    232211  inline G4bool GetApplyCutInRange(){ return ApplyCutInRange;}
    233   void  DefineCurrentMaterial(const G4MaterialCutsCouple* couple);
    234  
    235   inline G4String GetName(){ return name;}
    236 
    237 private: //Methods
    238  
    239  
    240  
    241                                    
    242 protected:
     212 
     213  inline G4String GetName(){ return name;}
     214  inline virtual void SetCSBiasingFactor(G4double aVal) {CS_biasing_factor = aVal;}
     215                                     
     216protected:
     217
     218  //Some of them can be overriden by daughter classes
     219 
     220 
     221  G4double DiffCrossSectionFunction1(G4double kinEnergyProj);
     222  G4double DiffCrossSectionFunction2(G4double kinEnergyProj);
     223  G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double EkinProd);
     224 
     225 
     226                               
     227  //General methods to sample secondary energy
     228  //--------------------------------------
     229  G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double prim_energy,G4bool IsScatProjToProjCase);
     230  G4double SampleAdjSecEnergyFromCSMatrix(G4double prim_energy,G4bool IsScatProjToProjCase);
     231  void     SelectCSMatrix(G4bool IsScatProjToProjCase);                       
     232 
     233  virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase);
     234 
     235 
     236 
     237  //Post  Step weight correction
     238  //----------------------------
     239  virtual void CorrectPostStepWeight(G4ParticleChange* fParticleChange,
     240                                     G4double old_weight,
     241                                     G4double adjointPrimKinEnergy,
     242                                     G4double projectileKinEnergy,
     243                                     G4bool IsScatProjToProjCase);     
     244 
     245 
     246 
     247 
     248 
     249 
     250protected: //attributes
    243251 
    244252  G4VEmModel* theDirectEMModel;
    … …  
    246254 
    247255
    248 protected:
    249 
    250   //  hide assignment operator
    251   G4VEmAdjointModel & operator=(const  G4VEmAdjointModel &right);
    252   G4VEmAdjointModel(const  G4VEmAdjointModel&);
    253  
    254  
     256
     257 
    255258  //Name
    256259  //-----
    … …  
    266269  G4double kinEnergyProdForIntegration;
    267270  G4double kinEnergyScatProjForIntegration;
     271  G4double kinEnergyProjForIntegration;
    268272 
    269273 
    … …  
    274278  std::vector< G4AdjointCSMatrix* >* pOnCSMatrixForProdToProjBackwardScattering;
    275279  std::vector< G4AdjointCSMatrix* >* pOnCSMatrixForScatProjToProjBackwardScattering;
    276   std::vector<double> CS_Vs_ElementForScatProjToProjCase;
    277   std::vector<double> CS_Vs_ElementForProdToProjCase;
     280  std::vector<G4double> CS_Vs_ElementForScatProjToProjCase;
     281  std::vector<G4double> CS_Vs_ElementForProdToProjCase;
    278282 
    279283  G4double lastCS;
     284  G4double lastAdjointCSForScatProjToProjCase;
     285  G4double lastAdjointCSForProdToProjCase;
     286 
     287 
     288 
    280289 
    281290  //particle definition
    … …  
    286295  G4ParticleDefinition* theDirectPrimaryPartDef;
    287296  G4bool second_part_of_same_type;
     297 
     298 
     299  //Prestep energy
     300  //-------------
     301  G4double preStepEnergy;
    288302 
    289303  //Current couple material
    … …  
    298312 
    299313 
    300   //CorrectWeightMode
    301   //------------------
    302  
    303   bool CorrectWeightMode;
    304  
    305   //Apply biasing
    306   //------------
    307  
    308   bool ApplyBiasing;
     314 
    309315 
    310316 
    … …  
    317323  G4double HighEnergyLimit;
    318324  G4double LowEnergyLimit;
     325
    319326 
    320327 
    … …  
    326333  //Type of Model with Matrix or not
    327334  //--------------------------------
    328    bool UseMatrix;
    329    bool UseMatrixPerElement; //other possibility is per Material
    330    bool UseOnlyOneMatrixForAllElements;
    331    bool IsIonisation;
     335   G4bool UseMatrix;
     336   G4bool UseMatrixPerElement; //other possibility is per Material
     337   G4bool UseOnlyOneMatrixForAllElements;
     338 
     339 
     340   //Index of Cross section matrices to be used
     341   //------------
     342   size_t indexOfUsedCrossSectionMatrix;
     343   
     344   
     345   
    332346   
    333347   

  • trunk/source/processes/electromagnetic/adjoint/include/G4eInverseBremsstrahlung.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseBremsstrahlung.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    2730//      Module:         G4eInverseBremstrahlung.hh
    2831//      Author:         L. Desorgher
    29 //      Date:           25 October 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    4649#define G4eInverseBremsstrahlung_h 1
    4750
    48 #include "G4VAdjointInverseScattering.hh"
     51#include "G4VAdjointReverseReaction.hh"
    4952#include "globals.hh"
    5053#include "G4eIonisation.hh"
    5154class G4AdjointBremsstrahlungModel;
    52 class G4eInverseBremsstrahlung: public G4VAdjointInverseScattering
     55class G4eInverseBremsstrahlung: public G4VAdjointReverseReaction
    5356
    5457{

  • trunk/source/processes/electromagnetic/adjoint/include/G4eInverseCompton.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseCompton.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    2730//      Module:         G4eInverseCompton.hh
    2831//      Author:         L. Desorgher
    29 //      Date:           25 October 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    4649#define G4eInverseCompton_h 1
    4750
    48 #include "G4VAdjointInverseScattering.hh"
     51#include "G4VAdjointReverseReaction.hh"
    4952#include "globals.hh"
    5053#include "G4eIonisation.hh"
    5154class G4AdjointComptonModel;
    52 class G4eInverseCompton: public G4VAdjointInverseScattering
     55class G4eInverseCompton: public G4VAdjointReverseReaction
    5356
    5457{

  • trunk/source/processes/electromagnetic/adjoint/include/G4eInverseIonisation.hh

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseIonisation.hh,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629/////////////////////////////////////////////////////////////////////////////////
    2730//      Module:         G4eInverseIonisation.hh
    2831//      Author:         L. Desorgher
    29 //      Date:           15 April 2007
    3032//      Organisation:   SpaceIT GmbH
     33//      Contract:       ESA contract 21435/08/NL/AT
    3134//      Customer:       ESA/ESTEC
    3235/////////////////////////////////////////////////////////////////////////////////
    … …  
    3942//-------------------------------------------------------------
    4043//      Documentation:
    41 //              Adjoint/revrese discrete ionisation
     44//              Adjoint/reverse discrete ionisation
    4245//
    4346
    … …  
    4548#define G4eInverseIonisation_h 1
    4649
    47 #include "G4VAdjointInverseScattering.hh"
     50#include "G4VAdjointReverseReaction.hh"
    4851#include "globals.hh"
    4952#include "G4eIonisation.hh"
    5053#include "G4VEmAdjointModel.hh"
    51 class G4eInverseIonisation: public G4VAdjointInverseScattering
     54class G4eInverseIonisation: public G4VAdjointReverseReaction
    5255
    5356{

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointAlongStepWeightCorrection.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointAlongStepWeightCorrection.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointAlongStepWeightCorrection.hh"
    2730#include "G4Step.hh"
    … …  
    3033#include "G4AdjointCSManager.hh"
    3134
     35#ifdef WIN32
     36#include <G4float .h>
     37#endif
    3238
    3339///////////////////////////////////////////////////////
    … …  
    4551{;
    4652}
    47 
    48 
    4953///////////////////////////////////////////////////////
    5054//
    51 
    5255void G4AdjointAlongStepWeightCorrection::PreparePhysicsTable(
    5356     const G4ParticleDefinition& )
    5457{
    5558;
    56 
    5759}
    58 
    5960///////////////////////////////////////////////////////
    6061//
    … …  
    6364{;
    6465}
    65 
    66 
    67 
    68 
    6966///////////////////////////////////////////////////////
    7067//
    … …  
    8582                                                                        preStepKinEnergy,Tkin, currentCouple,length);
    8683       
     84 
    8785 
     86  G4double new_weight=weight_correction*track.GetWeight();
    8887 
    89   G4double new_weight=weight_correction*track.GetWeight();
     88  //if (weight_correction >2.) new_weight=1.e-300;
     89 
     90 
     91  //The following test check for zero weight.
     92  //This happens after weight correction of gamma for photo electric effect.
     93  //When the new weight is 0 it will be later on consider as nan by G4.
     94  //Therefore we do put a lower limit of 1.e-300. for new_weight
     95  //Correction by L.Desorgher on 15 July 2009
     96#ifdef WIN32
     97  if (!!_isnan(new_weight) || new_weight==0){
     98#else
     99  if (std::isnan(new_weight) || new_weight==0){
     100#endif
     101                //G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
     102                new_weight=1.e-300;
     103  }
     104 
     105  //G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
    90106  fParticleChange->SetParentWeightByProcess(false);
    91107  fParticleChange->SetSecondaryWeightByProcess(false);
    … …  
    96112
    97113}
     114///////////////////////////////////////////////////////
     115//
     116G4double G4AdjointAlongStepWeightCorrection::GetContinuousStepLimit(const G4Track& track,
     117                G4double , G4double , G4double& )
     118{
     119  G4double x = DBL_MAX;
     120  DefineMaterial(track.GetMaterialCutsCouple());
     121  preStepKinEnergy = track.GetKineticEnergy();
     122  return x;
     123}

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointBremsstrahlungModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointBremsstrahlungModel.hh"
    2730#include "G4AdjointCSManager.hh"
    … …  
    3033#include "G4ParticleChange.hh"
    3134#include "G4AdjointElectron.hh"
     35#include "G4AdjointGamma.hh"
     36#include "G4Electron.hh"
     37
    3238#include "G4Timer.hh"
     39//#include "G4PenelopeBremsstrahlungModel.hh"
     40
    3341
    3442////////////////////////////////////////////////////////////////////////////////
    3543//
    3644G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
    37  G4VEmAdjointModel("AdjointBremModel"),
    38  probsup(1.0),
    39  MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
    40  LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
    41  theLPMflag(true)
    42 
    43 { isElectron= true;
    44   SetUseMatrix(true);
     45 G4VEmAdjointModel("AdjointeBremModel"),
     46  MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi)
     47{
     48  SetUseMatrix(false);
    4549  SetUseMatrixPerElement(false);
     50 
     51  theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel");
     52  theDirectEMModel=theDirectStdBremModel;
     53 // theDirectPenelopeBremModel =0;
     54       
    4655  SetApplyCutInRange(true);
    47   SetIsIonisation(false);
    4856  highKinEnergy= 100.*TeV;
    4957  lowKinEnergy = 1.0*keV;
    5058  theTimer =new G4Timer();
    5159 
    52   theTimer->Start();
    53   InitialiseParameters();
    54   theTimer->Stop();
    55   G4cout<<"Time elapsed in second for the initialidation of AdjointBrem "<<theTimer->GetRealElapsed()<<std::endl;
    56  
    57   ModeldCS="MODEL1";
     60  theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
     61  theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
     62  theDirectPrimaryPartDef=G4Electron::Electron();
     63  second_part_of_same_type=false;
     64 
     65  /*UsePenelopeModel=false;
     66  if (UsePenelopeModel) {
     67        G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
     68        theEmModelManagerForFwdModels = new G4EmModelManager();
     69        isPenelopeModelInitialised = false;
     70        G4VEmFluctuationModel* f=0;
     71        G4Region* r=0;
     72        theDirectEMModel=thePenelopeModel;
     73        theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
     74  }
     75  */   
     76 
     77
    5878 
    5979}
     80
     81////////////////////////////////////////////////////////////////////////////////
     82//
     83void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
     84                       G4bool IsScatProjToProjCase,
     85                       G4ParticleChange* fParticleChange)
     86{
     87 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
     88
     89 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
     90 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
     91 
     92 
     93 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
     94 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
     95 
     96 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
     97        return;
     98 }
     99 
     100  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
     101                                                        IsScatProjToProjCase);
     102 //Weight correction
     103 //-----------------------                                         
     104 CorrectPostStepWeight(fParticleChange,
     105                       aTrack.GetWeight(),
     106                       adjointPrimKinEnergy,
     107                       projectileKinEnergy,
     108                       IsScatProjToProjCase);   
     109 
     110 
     111 //Kinematic
     112 //---------
     113 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
     114 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
     115 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
     116 G4double projectileP = std::sqrt(projectileP2);
     117 
     118 
     119 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
     120 //------------------------------------------------
     121  G4double u;
     122  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
     123
     124  if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
     125     else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
     126
     127  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
     128
     129  G4double sint = std::sin(theta);
     130  G4double cost = std::cos(theta);
     131
     132  G4double phi = twopi * G4UniformRand() ;
     133 
     134  G4ThreeVector projectileMomentum;
     135  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
     136  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
     137        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
     138        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
     139        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
     140        G4double sint1 =  std::sqrt(1.-cost1*cost1);
     141        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
     142 
     143  }
     144 
     145  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
     146 
     147 
     148 
     149  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
     150        fParticleChange->ProposeTrackStatus(fStopAndKill);
     151        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
     152  }
     153  else {
     154        fParticleChange->ProposeEnergy(projectileKinEnergy);
     155        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
     156       
     157  }     
     158}
     159////////////////////////////////////////////////////////////////////////////////
     160//
     161void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
     162                       G4bool IsScatProjToProjCase,
     163                       G4ParticleChange* fParticleChange)
     164{
     165
     166 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
     167 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
     168 
     169 
     170 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
     171 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
     172 
     173 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
     174        return;
     175 }
     176 
     177 G4double projectileKinEnergy =0.;
     178 G4double gammaEnergy=0.;
     179 G4double diffCSUsed=0.;
     180 if (!IsScatProjToProjCase){
     181        gammaEnergy=adjointPrimKinEnergy;
     182        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
     183        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
     184        if (Emin>=Emax) return;
     185        projectileKinEnergy=Emin*pow(Emax/Emin,G4UniformRand());
     186        diffCSUsed=lastCZ/projectileKinEnergy;
     187       
     188 }
     189 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
     190        G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
     191        if (Emin>=Emax) return;
     192        G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
     193        G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
     194        //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
     195        projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*pow(f2,G4UniformRand()));
     196        gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
     197        diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
     198       
     199 }
     200 
     201 
     202 
     203                                                       
     204 //Weight correction
     205 //-----------------------
     206 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
     207 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
     208
     209 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
     210 //Here we consider the true  diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term.
     211 //Basically any other differential CS   diffCS could be used here (example Penelope).
     212 
     213 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
     214 w_corr*=diffCS/diffCSUsed;
     215           
     216 G4double new_weight = aTrack.GetWeight()*w_corr;
     217 fParticleChange->SetParentWeightByProcess(false);
     218 fParticleChange->SetSecondaryWeightByProcess(false);
     219 fParticleChange->ProposeParentWeight(new_weight);
     220 
     221 //Kinematic
     222 //---------
     223 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
     224 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
     225 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
     226 G4double projectileP = std::sqrt(projectileP2);
     227 
     228 
     229 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
     230 //------------------------------------------------
     231  G4double u;
     232  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
     233
     234  if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
     235     else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
     236
     237  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
     238
     239  G4double sint = std::sin(theta);
     240  G4double cost = std::cos(theta);
     241
     242  G4double phi = twopi * G4UniformRand() ;
     243 
     244  G4ThreeVector projectileMomentum;
     245  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
     246  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
     247        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
     248        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
     249        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
     250        G4double sint1 =  std::sqrt(1.-cost1*cost1);
     251        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
     252 
     253  }
     254 
     255  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
     256 
     257 
     258 
     259  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
     260        fParticleChange->ProposeTrackStatus(fStopAndKill);
     261        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
     262  }
     263  else {
     264        fParticleChange->ProposeEnergy(projectileKinEnergy);
     265        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
     266       
     267  }     
     268}
    60269////////////////////////////////////////////////////////////////////////////////
    61270//
    62271G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
    63272{;}
    64 ////////////////////////////////////////////////////////////////////////////////
    65 //
    66 /*G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(
    67                                       const G4Material* aMaterial,
    68                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    69                                       G4double kinEnergyProd // kinetic energy of the secondary particle
    70                                       )
    71 {
    72  
    73  static const G4double
    74      ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
    75      ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
    76      ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
    77 
    78   static const G4double
    79      bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
    80      bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
    81      bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
    82 
    83   static const G4double
    84      al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
    85      al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
    86      al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
    87 
    88   static const G4double
    89      bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
    90      bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
    91      bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
    92 
    93   static const G4double tlow = 1.*MeV;
    94  
    95   G4double dCrossEprod=0.;
    96   G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
    97   G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
    98  
    99  
    100  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
    101        
    102         G4double cross = 0.0;
    103        
    104 
    105        
    106         G4double E1=kinEnergyProd;
    107         G4double E2=kinEnergyProd*1.000000001;
    108         G4double dE=(E2-E1);
    109 
    110         const G4ElementVector* theElementVector = aMaterial->GetElementVector();
    111         const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
    112         G4double dum=0.;
    113  
    114         for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
    115 
    116                 G4double fac=
    117                
    118                 cross += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
    119                       kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
    120        
    121                
    122                
    123         }
    124         dCrossEprod=(cross1-cross2)/dE; //first term
    125        
    126         //Now come the correction
    127         //-----------------------
    128        
    129         //First compute fsig for E1
    130         //-------------------------
    131        
    132        
    133         G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
    134         G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
    135                                              *(aMaterial->GetElectronDensity());
    136 
    137         G4double fsig = 0.;
    138         G4int nmax = 100;
    139         G4double vmin=std::log(E1);
    140         G4double vmax=std::log(kinEnergyProj) ;
    141         G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
    142         G4double u,fac,c,v,dv,y ;
    143         if(nn > 0) {
    144 
    145                 dv = (vmax-vmin)/nn ;
    146                 v  = vmin-dv ;
    147                 for(G4int n=0; n<=nn; n++) {
    148 
    149                         v += dv; 
    150                         u = std::exp(v);             
    151                         fac = SupressionFunction(aMaterial, kinEnergyProj, u);
    152                         y = u/kinEnergyProj;
    153                         fac *= (4.-4.*y+3.*y*y)/3.;
    154                         fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
    155 
    156                         if ((n==0)||(n==nn)) c=0.5;
    157                         else    c=1. ;
    158 
    159                         fac  *= c;
    160                         fsig += fac;
    161                 }
    162                 y = E1/kinEnergyProj ;
    163                 fsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
    164 
    165         }
    166         else {
    167                 fsig = 1.;
    168         }
    169         if (fsig > 1.) fsig = 1.;
    170        
    171         dCrossEprod*=fsig;
    172         //return dCrossEprod;
    173         //Now we  compute dfsig
    174         //-------------------------
    175         G4double dfsig = 0.;
    176         nn=20;
    177         vmax=std::log(E2) ;
    178         dv = (vmax-vmin)/nn ;
    179         v  = vmin-dv ;
    180         for(G4int n=0; n<=nn; n++) {
    181                 v += dv; 
    182                 u = std::exp(v);             
    183                 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
    184                 y = u/kinEnergyProj;
    185                 fac *= (4.-4.*y+3.*y*y)/3.;
    186                 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
    187 
    188                 if ((n==0)||(n==nn)) c=0.5;
    189                 else    c=1. ;
    190 
    191                 fac  *= c;
    192                 dfsig += fac;
    193         }
    194         y = E1/kinEnergyProj;
    195         dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
    196        
    197         dCrossEprod+=dfsig*cross1/dE;
    198        
    199        
    200        
    201          
    202        
    203  }
    204  return dCrossEprod;
    205  
    206 }
    207 */
     273
     274////////////////////////////////////////////////////////////////////////////////
     275//
    208276G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
    209277                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    210278                                      G4double kinEnergyProd // kinetic energy of the secondary particle
    211279                                      )
    212 {if (ModeldCS=="MODEL2") return  DiffCrossSectionPerVolumePrimToSecond2(aMaterial,
    213                                                                  kinEnergyProj,  // kinetic energy of the primary particle before the interaction
     280{/*if (UsePenelopeModel && !isPenelopeModelInitialised) {
     281        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
     282        isPenelopeModelInitialised =true;
     283 }
     284 */                                                             
     285 return  DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
     286                                                                 kinEnergyProj,
    214287                                                                 kinEnergyProd);
    215  if (ModeldCS=="MODEL3") return  DiffCrossSectionPerVolumePrimToSecond3(aMaterial,
    216                                                                  kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    217                                                                  kinEnergyProd);
    218  return  DiffCrossSectionPerVolumePrimToSecond1(aMaterial,
    219                                                                  kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    220                                                                  kinEnergyProd);                                                                 
     288 /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
     289                                                                 kinEnergyProj,
     290                                                                 kinEnergyProd);*/                                                                                                                               
    221291}                                     
    222 ////////////////////////////////////////////////////////////////////////////////
    223 // the one used till now
    224 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond1(
     292
     293////////////////////////////////////////////////////////////////////////////////
     294//
     295G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
    225296                                      const G4Material* aMaterial,
    226297                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    … …  
    233304 
    234305 
     306 //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
     307 //This is what is applied in the discrete standard model before the  rejection test  that make a cooerction
     308 //The application of the same rejection function is not possble here.
     309 //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the
     310 // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and
     311 // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one.
     312 // In the future we plan to use the brem secondary spectra from the G4Penelope implementation
     313 
    235314 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
    236        
    237         G4double cross1 = 0.0;
    238         G4double cross2 = 0.0;
    239 
    240        
    241         G4double E1=kinEnergyProd;
    242         G4double E2=kinEnergyProd*1.01;
    243         G4double dE=(E2-E1);
    244 
    245         const G4ElementVector* theElementVector = aMaterial->GetElementVector();
    246         const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
    247         G4double dum=0.;
    248  
    249         for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
    250 
    251                 cross1 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
    252                       kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
    253        
    254                 cross2 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
    255                      kinEnergyProj, (*theElementVector)[i]->GetZ(), dum, E2);
    256                
    257         }
    258         dCrossEprod=(cross1-cross2)/dE; //first term
    259        
    260         //Now come the correction
    261         //-----------------------
    262        
    263         //First compute fsig for E1
    264         //-------------------------
    265        
    266        
    267         G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
    268         G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
    269                                              *(aMaterial->GetElectronDensity());
    270 
    271         G4double fsig1 = 0.;
    272         G4int nmax = 100;
    273         G4double vmin=std::log(E1);
    274         G4double vmax=std::log(kinEnergyProj) ;
    275         G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
    276         G4double u,fac,c,v,dv,y ;
    277         if(nn > 0) {
    278 
    279                 dv = (vmax-vmin)/nn ;
    280                 v  = vmin-dv ;
    281                 for(G4int n=0; n<=nn; n++) {
    282 
    283                         v += dv; 
    284                         u = std::exp(v);             
    285                         fac = SupressionFunction(aMaterial, kinEnergyProj, u);
    286                         y = u/kinEnergyProj;
    287                         fac *= (4.-4.*y+3.*y*y)/3.;
    288                         fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
    289 
    290                         if ((n==0)||(n==nn)) c=0.5;
    291                         else    c=1. ;
    292 
    293                         fac  *= c;
    294                         fsig1 += fac;
    295                 }
    296                 y = E1/kinEnergyProj ;
    297                 fsig1 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
    298 
    299         }
    300         else {
    301                 fsig1 = 1.;
    302         }
    303         if (fsig1 > 1.) fsig1 = 1.;
    304        
    305         dCrossEprod*=fsig1;
    306        
    307        
    308         G4double fsig2 = 0.;
    309         vmin=std::log(E2);
    310         nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
    311         if(nn > 0) {
    312 
    313                 dv = (vmax-vmin)/nn ;
    314                 v  = vmin-dv ;
    315                 for(G4int n=0; n<=nn; n++) {
    316 
    317                         v += dv; 
    318                         u = std::exp(v);             
    319                         fac = SupressionFunction(aMaterial, kinEnergyProj, u);
    320                         y = u/kinEnergyProj;
    321                         fac *= (4.-4.*y+3.*y*y)/3.;
    322                         fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
    323 
    324                         if ((n==0)||(n==nn)) c=0.5;
    325                         else    c=1. ;
    326 
    327                         fac  *= c;
    328                         fsig2 += fac;
    329                 }
    330                 y = E2/kinEnergyProj ;
    331                 fsig2 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
    332 
    333         }
    334         else {
    335                 fsig2 = 1.;
    336         }
    337         if (fsig2 > 1.) fsig2 = 1.;
    338        
    339 
    340         G4double dfsig=(fsig2-fsig1);
    341         dCrossEprod+=dfsig*cross1/dE;
    342        
    343         dCrossEprod=(fsig1*cross1-fsig2*cross2)/dE;
    344        
    345        
    346        
    347        
    348        
    349         /*if (fsig < 1.){
    350                 //Now we  compute dfsig
    351                 //-------------------------
    352                 G4double dfsig = 0.;
    353                 nn=20;
    354                 vmax=std::log(E2) ;
    355                 dv = (vmax-vmin)/nn ;
    356                 v  = vmin-dv ;
    357                 for(G4int n=0; n<=nn; n++) {
    358                         v += dv; 
    359                         u = std::exp(v);             
    360                         fac = SupressionFunction(aMaterial, kinEnergyProj, u);
    361                         y = u/kinEnergyProj;
    362                         fac *= (4.-4.*y+3.*y*y)/3.;
    363                         fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
    364 
    365                         if ((n==0)||(n==nn)) c=0.5;
    366                         else    c=1. ;
    367 
    368                         fac  *= c;
    369                         dfsig += fac;
    370                 }
    371                 y = E1/kinEnergyProj;
    372                 dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
    373                 dCrossEprod+=dfsig*cross1/dE;
    374                
    375         }       
    376         */
    377        
    378        
    379        
    380        
    381        
    382          
    383        
     315        G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
     316        dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
    384317 }
    385318 return dCrossEprod;
    386319 
    387 }
    388 
    389 
    390 ////////////////////////////////////////////////////////////////////////////////
    391 //
    392 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2(
    393                                       const G4Material* aMaterial,
     320}
     321
     322////////////////////////////////////////////////////////////////////////////////
     323//
     324G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
     325                                      const G4Material* material,
    394326                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    395327                                      G4double kinEnergyProd // kinetic energy of the secondary particle
    396328                                      )
    397329{
     330 //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor
     331  //used in the direct model
     332 
    398333 G4double dCrossEprod=0.;
    399  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
    400  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
    401  
    402  
    403  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
    404        
    405         G4double dEdX1 = 0.0;
    406         G4double dEdX2 = 0.0;
    407 
    408        
    409         G4double E1=kinEnergyProd;
    410         G4double E2=kinEnergyProd*1.001;
    411         G4double dE=(E2-E1);
    412         //G4double dum=0.;
    413        
    414         dEdX1 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E1);
    415         dEdX2 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E2);
    416         dCrossEprod=(dEdX2-dEdX1)/dE/E1;
    417        
    418          
    419        
    420        
    421        
    422          
    423        
    424  }
     334 
     335 const G4ElementVector* theElementVector = material->GetElementVector();
     336 const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
     337 G4double dum=0.;
     338 G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
     339 G4double dE=E2-E1;
     340 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
     341        G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
     342        G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
     343        dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
     344   
     345 }
     346 
     347 //Now the Migdal correction
     348 
     349 G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
     350 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
     351                                             *(material->GetElectronDensity());
     352                                             
     353 
     354 G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
     355                                                                    //model is different than the one used in the secondary sampling by a
     356                                                                    //factor (1.+kp2) To be checked!
     357 
     358 dCrossEprod*=MigdalFactor;
    425359 return dCrossEprod;
    426360 
    … …  
    428362////////////////////////////////////////////////////////////////////////////////
    429363//
    430 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond3(
    431                                       const G4Material* aMaterial,
    432                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    433                                       G4double kinEnergyProd // kinetic energy of the secondary particle
    434                                       )
    435 {
    436  
    437  return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
    438                                                                  kinEnergyProj,  // kinetic energy of the primary particle before the interaction
    439                                                                  kinEnergyProd);
    440  
    441 }
    442 
    443 ////////////////////////////////////////////////////////////////////////////////
    444 //
    445 G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material,
    446                                  G4double kineticEnergy, G4double gammaEnergy)
    447 {
    448   // supression due to the LPM effect+polarisation of the medium/
    449   // supression due to the polarisation alone
    450 
    451 
    452   G4double totEnergy = kineticEnergy+electron_mass_c2 ;
    453   G4double totEnergySquare = totEnergy*totEnergy ;
    454 
    455   G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
    456 
    457   G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
    458 
    459   G4double electronDensity = material->GetElectronDensity();
    460 
    461   G4double sp = gammaEnergySquare/
    462    (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
    463 
    464   G4double supr = 1.0;
    465 
    466   if (theLPMflag) {
    467 
    468     G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
    469 
    470     if (s2lpm < 1.) {
    471 
    472       G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
    473       G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
    474       G4double splim = LPMgEnergyLimit2/
    475         (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
    476       G4double w = 1.+1./splim ;
    477 
    478       if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
    479       else                 w = s2lpm*(1.+1./sp);
    480 
    481       supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ;
    482       supr /= sp;   
    483     }
    484    
    485   }
    486   return supr;
    487 }
    488 
    489 ////////////////////////////////////////////////////////////////////////////////
    490 //
    491 void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
    492                        G4bool IsScatProjToProjCase,
    493                        G4ParticleChange* fParticleChange)
    494 {
    495 
    496   //G4cout<<"Adjoint Brem"<<std::endl;
    497   const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
    498  
    499   size_t ind=0;
    500  
    501   if (UseMatrixPerElement ) { //Select Material
    502         std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
    503         if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
    504         G4double rand_var= G4UniformRand();
    505         G4double SumCS=0.;
    506         for (size_t i=0;i<CS_Vs_Element->size();i++){
    507                 SumCS+=(*CS_Vs_Element)[i];
    508                 if (rand_var<=SumCS/lastCS){
    509                         ind=i;
    510                         break;
    511                 }
    512         }
    513   }
    514   else  {
    515         ind = currentMaterialIndex;
    516   }
    517  
    518  
    519  //Elastic inverse scattering modified compared to general G4VEmAdjointModel
    520  //---------------------------
    521  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
    522  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
    523  //G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
    524  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
    525         return;
    526  }
    527  
    528  //Sample secondary energy
    529  //-----------------------
    530  
    531  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
    532                                                    adjointPrimKinEnergy,
    533                                                    IsScatProjToProjCase);
    534                                    
    535  
    536  
    537  
    538  //Weight correction
    539  //-----------------------                                         
    540  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy); 
    541  
    542  
    543  //Kinematic
    544  //---------
    545  
    546  G4double projectileM0 = electron_mass_c2;
    547  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
    548  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
    549  G4double projectileP = std::sqrt(projectileP2);
    550  
    551  
    552  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
    553  //------------------------------------------------
    554   G4double u;
    555   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
    556 
    557   if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
    558      else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
    559 
    560   G4double theta = u*electron_mass_c2/projectileTotalEnergy;
    561 
    562   G4double sint = std::sin(theta);
    563   G4double cost = std::cos(theta);
    564 
    565   G4double phi = twopi * G4UniformRand() ;
    566  
    567   G4ThreeVector projectileMomentum;
    568   projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
    569   if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
    570         G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
    571         G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
    572         G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
    573         G4double sint1 =  std::sqrt(1.-cost1*cost1);
    574         projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
    575  
    576   }
    577  
    578   projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
    579  
    580  
    581  
    582   if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
    583         fParticleChange->ProposeTrackStatus(fStopAndKill);
    584         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
    585         //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
     364G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
     365                                             G4double primEnergy,
     366                                             G4bool IsScatProjToProjCase)
     367{/* if (UsePenelopeModel && !isPenelopeModelInitialised) {
     368        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
     369        isPenelopeModelInitialised =true;
     370  }
     371  */
     372  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
     373  DefineCurrentMaterial(aCouple);
     374  G4double Cross=0.;
     375  lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/exp(1.));//this give the constant above
     376 
     377  if (!IsScatProjToProjCase ){
     378        G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
     379        G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
     380        if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*log(Emax_proj/Emin_proj);
    586381  }
    587382  else {
    588         fParticleChange->ProposeEnergy(projectileKinEnergy);
    589         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
    590         //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
    591   }     
    592 }
    593 ////////////////////////////////////////////////////////////////////////////////
    594 //
    595 void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel)
    596 {theDirectBremModel=aModel;
    597  DefineDirectEMModel(aModel);
    598 }
    599 ////////////////////////////////////////////////////////////////////////////////
    600 //
    601 void G4AdjointBremsstrahlungModel::InitialiseParameters()
    602 { 
    603   static const G4double
    604      ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
    605      ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
    606      ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
    607 
    608   static const G4double
    609      bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
    610      bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
    611      bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
    612 
    613  /* static const G4double
    614      al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
    615      al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
    616      al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
    617 
    618   static const G4double
    619      bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
    620      bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
    621      bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/
    622  
    623  
    624   const G4ElementTable* theElementTable = G4Element::GetElementTable();
    625   FZ.clear();
    626   ah1.clear();
    627   ah2.clear();
    628   ah3.clear();
    629  
    630   bh1.clear();
    631   bh2.clear();
    632   bh3.clear();
    633  
    634   al0.clear();
    635   al1.clear();
    636   al2.clear();
    637  
    638   bl0.clear();
    639   bl1.clear();
    640   bl2.clear();
    641   SigmaPerAtom.clear();
    642  
    643   for (size_t j=0; j<theElementTable->size();j++){
    644        
    645         G4Element* anElement=(*theElementTable)[j];
    646         G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
    647         FZ.push_back(lnZ* (4.- 0.55*lnZ));
    648         G4double ZZ = anElement->GetIonisation()->GetZZ3();
    649        
    650         ah1.push_back(ah10 + ZZ* (ah11 + ZZ* ah12));
    651         ah2.push_back(ah20 + ZZ* (ah21 + ZZ* ah22));
    652         ah3.push_back(ah30 + ZZ* (ah31 + ZZ* ah32));
    653 
    654         bh1.push_back(bh10 + ZZ* (bh11 + ZZ* bh12));
    655         bh2.push_back(bh20 + ZZ* (bh21 + ZZ* bh22));
    656         bh3.push_back(bh30 + ZZ* (bh31 + ZZ* bh32));
    657         /*SigmaPerAtom.push_back(theDirectEMModel->ComputeCrossSectionPerAtom(
    658                                         theDirectPrimaryPartDef,GetHighEnergyLimit()/2.,
    659                                         anElement->GetZ(),1.,GetLowEnergyLimit(),1.e20));*/
    660        
    661        
     383        G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
     384        G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
     385        if (Emax_proj>Emin_proj) Cross= lastCZ*log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
    662386       
    663   }     
    664 }
     387  }
     388  return Cross;
     389}                                           
     390
     391
     392
     393
     394

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointCSManager.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointCSManager.hh"
    2730#include "G4AdjointCSMatrix.hh"
    … …  
    4043#include "G4Electron.hh"
    4144#include "G4Gamma.hh"
     45#include "G4Proton.hh"
    4246#include "G4AdjointElectron.hh"
    4347#include "G4AdjointGamma.hh"
     48#include "G4AdjointProton.hh"
    4449#include "G4ProductionCutsTable.hh"
    4550#include "G4ProductionCutsTable.hh"
     51#include <fstream>
     52#include <iomanip>
    4653
    4754
    … …  
    6572  listOfForwardEmProcess.clear();
    6673  listOfForwardEnergyLossProcess.clear();
    67   theListOfAdjointParticlesInAction.clear();
     74  theListOfAdjointParticlesInAction.clear();
     75  EminForFwdSigmaTables.clear();
     76  EminForAdjSigmaTables.clear();
     77  EkinofFwdSigmaMax.clear();
     78  EkinofAdjSigmaMax.clear();
    6879  Tmin=0.1*keV;
    6980  Tmax=100.*TeV;
    70   nbins=240;
     81  nbins=360; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
    7182 
    7283  RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
    7384  RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
     85  RegisterAdjointParticle(G4AdjointProton::AdjointProton());
    7486 
    7587  verbose  = 1;
    76  
    77   consider_continuous_weight_correction =true;
    78   consider_poststep_weight_correction =false;
     88 
     89  lastPartDefForCS =0;
     90  LastEkinForCS =0;
     91  LastCSCorrectionFactor =1.;
     92 
     93  forward_CS_mode = true;
     94 
     95  currentParticleDef = 0;
     96 
     97  theAdjIon = 0;
     98  theFwdIon = 0; 
     99 
    79100 
    80101}
    … …  
    96117  if (anAdjPartDef && aProcess){
    97118        RegisterAdjointParticle(anAdjPartDef);
    98         int index=-1;
     119        G4int index=-1;
    99120       
    100121        for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    … …  
    111132  if (anAdjPartDef && aProcess){
    112133        RegisterAdjointParticle(anAdjPartDef);
    113         int index=-1;
     134        G4int index=-1;
    114135        for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    115136                if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
    … …  
    121142//
    122143void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
    123 {  int index=-1;
     144{  G4int index=-1;
    124145   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    125146        if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
    … …  
    132153        listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
    133154        theListOfAdjointParticlesInAction.push_back(aPartDef);
     155        EminForFwdSigmaTables.push_back(std::vector<G4double> ());
     156        EminForAdjSigmaTables.push_back(std::vector<G4double> ());
     157        EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
     158        EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
     159       
    134160   }
    135161}
    … …  
    162188        const G4ElementTable* theElementTable = G4Element::GetElementTable();
    163189        const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
     190       
     191        G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
    164192        for (size_t i=0; i<listOfAdjointEMModel.size();i++){
    165193                G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
    166                 G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<std::endl;
     194                G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
    167195                if (aModel->GetUseMatrix()){
    168196                        std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
    … …  
    173201                                if (aModel->GetUseOnlyOneMatrixForAllElements()){
    174202                                                std::vector<G4AdjointCSMatrix*>
    175                                                 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 10);
     203                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
    176204                                                aListOfMat1->push_back(two_matrices[0]);
    177205                                                aListOfMat2->push_back(two_matrices[1]);
    … …  
    180208                                        for (size_t j=0; j<theElementTable->size();j++){
    181209                                                G4Element* anElement=(*theElementTable)[j];
    182                                                 G4int Z = G4int(anElement->GetZ());
    183                                                 G4int A = G4int(anElement->GetA());
     210                                                G4int Z = int(anElement->GetZ());
     211                                                G4int A = int(anElement->GetA());
    184212                                                std::vector<G4AdjointCSMatrix*>
    185                                                         two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 10);
     213                                                        two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
    186214                                                aListOfMat1->push_back(two_matrices[0]);
    187215                                                aListOfMat2->push_back(two_matrices[1]);
    … …  
    193221                                        G4Material* aMaterial=(*theMaterialTable)[j];
    194222                                        std::vector<G4AdjointCSMatrix*>
    195                                                 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 10);
     223                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
    196224                                        aListOfMat1->push_back(two_matrices[0]);
    197225                                        aListOfMat2->push_back(two_matrices[1]);
    … …  
    203231                        aModel->SetCSMatrices(aListOfMat1, aListOfMat2);       
    204232                }
    205                 else {  std::vector<G4AdjointCSMatrix*> two_empty_matrices;
     233                else {  G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
     234                        std::vector<G4AdjointCSMatrix*> two_empty_matrices;
    206235                        theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
    207236                        theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
    … …  
    209238                }               
    210239        }
    211         G4cout<<"All adjoint cross section matrices are built "<<std::endl;
     240        G4cout<<"              All adjoint cross section matrices are computed!"<<G4endl;
     241        G4cout<<"======================================================================"<<G4endl;
     242       
    212243        CrossSectionMatrixesAreBuilt = true;
     244
     245
    213246}
    214247
    … …  
    221254  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    222255        G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
     256        DefineCurrentParticle(thePartDef);
    223257        theTotalForwardSigmaTableVector[i]->clearAndDestroy();
    224258        theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
     259        EminForFwdSigmaTables[i].clear();
     260        EminForAdjSigmaTables[i].clear();
     261        EkinofFwdSigmaMax[i].clear();
     262        EkinofAdjSigmaMax[i].clear();
     263        //G4cout<<thePartDef->GetParticleName();
     264       
    225265        for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
    226266                const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
    227267               
     268                /*
     269                G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
     270                G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
     271               
     272                std::fstream FileOutputAdjCS(file_name1, std::ios::out);
     273                std::fstream FileOutputFwdCS(file_name2, std::ios::out);
     274               
     275               
     276               
     277                FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
     278                FileOutputAdjCS<<std::setprecision(6);
     279                FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
     280                FileOutputFwdCS<<std::setprecision(6); 
     281                */
     282                         
     283
    228284                //make first the total fwd CS table for FwdProcess
    229285                G4PhysicsVector* aVector =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
     286                G4bool Emin_found=false;
     287                size_t ind=0;
     288                G4double sigma_max =0.;
     289                G4double e_sigma_max =0.;
    230290                for(size_t l=0; l<aVector->GetVectorLength(); l++) {
    231                         G4double totCS=0;
     291                        G4double totCS=0.;
    232292                        G4double e=aVector->GetLowEdgeEnergy(l);
    233293                        for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
    … …  
    235295                        }
    236296                        for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
    237                                 totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e, couple);
    238                         }
    239                         //G4cout<<totCS<<std::endl;
     297                                if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
     298                                        size_t mat_index = couple->GetIndex();
     299                                        G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
     300                                        G4double chargeSqRatio =  currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
     301                                        (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
     302                                }
     303                                G4double e1=e/massRatio;
     304                                totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
     305                        }
    240306                        aVector->PutValue(l,totCS);
    241                
    242                 }
     307                        if (totCS>sigma_max){
     308                                sigma_max=totCS;
     309                                e_sigma_max = e;
     310                               
     311                        }
     312                        //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
     313                       
     314                        if (totCS>0 && !Emin_found) {
     315                                EminForFwdSigmaTables[i].push_back(e);
     316                                Emin_found=true;
     317                        }
     318                       
     319               
     320                }
     321                //FileOutputFwdCS.close();
     322               
     323                EkinofFwdSigmaMax[i].push_back(e_sigma_max);
     324               
     325               
     326                if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
     327               
    243328                theTotalForwardSigmaTableVector[i]->push_back(aVector);
    244329               
     330               
     331                Emin_found=false;
     332                sigma_max=0;
     333                e_sigma_max =0.;
     334                ind=0;
    245335                G4PhysicsVector* aVector1 =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
    246336                for(size_t l=0; l<aVector->GetVectorLength(); l++) {
    247337                        G4double e=aVector->GetLowEdgeEnergy(l);
    248                         G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e);
    249                         //G4cout<<totCS<<std::endl;
     338                        G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
    250339                        aVector1->PutValue(l,totCS);
    251                
    252                 }       
     340                        if (totCS>sigma_max){
     341                                sigma_max=totCS;
     342                                e_sigma_max = e;
     343                               
     344                        }
     345                        //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
     346                        if (totCS>0 && !Emin_found) {
     347                                EminForAdjSigmaTables[i].push_back(e);
     348                                Emin_found=true;
     349                        }
     350               
     351                }
     352                //FileOutputAdjCS.close();
     353                EkinofAdjSigmaMax[i].push_back(e_sigma_max);
     354                if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
     355                       
    253356                theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
    254357               
    … …  
    262365                                     const G4MaterialCutsCouple* aCouple)
    263366{ DefineCurrentMaterial(aCouple);
    264   int index=-1;
    265   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    266         if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
    267   }     
    268   if (index == -1) return 0.;
    269  
     367  DefineCurrentParticle(aPartDef);     
    270368  G4bool b;
    271   return (((*theTotalAdjointSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
     369  return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
    272370 
    273371 
    … …  
    279377                                     const G4MaterialCutsCouple* aCouple)
    280378{ DefineCurrentMaterial(aCouple);
    281   int index=-1;
    282   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
    283         if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
    284   }     
    285   if (index == -1) return 0.;
     379  DefineCurrentParticle(aPartDef);
    286380  G4bool b;
    287   return (((*theTotalForwardSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
    288  
    289  
     381  return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
     382 
     383 
     384}
     385                                     
     386///////////////////////////////////////////////////////
     387//                                   
     388void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
     389                                     const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
     390{ DefineCurrentMaterial(aCouple);
     391  DefineCurrentParticle(aPartDef);
     392  emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
     393  emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
     394 
     395 
     396 
     397}
     398///////////////////////////////////////////////////////
     399//                                   
     400void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
     401                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
     402{ DefineCurrentMaterial(aCouple);
     403  DefineCurrentParticle(aPartDef);
     404  e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
     405  G4bool b;
     406  sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
     407  e_sigma_max/=massRatio;
     408 
     409 
     410}
     411///////////////////////////////////////////////////////
     412//                                   
     413void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
     414                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
     415{ DefineCurrentMaterial(aCouple);
     416  DefineCurrentParticle(aPartDef);
     417  e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
     418  G4bool b;
     419  sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
     420  e_sigma_max/=massRatio;
     421 
     422 
     423}                                   
     424///////////////////////////////////////////////////////
     425//                                   
     426G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
     427                                                                                G4double& fwd_TotCS)
     428{ G4double corr_fac = 1.;
     429  if (forward_CS_mode) {
     430        fwd_TotCS=PrefwdCS;
     431        if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
     432                DefineCurrentMaterial(aCouple);
     433                PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
     434                PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
     435                LastEkinForCS = PreStepEkin;
     436                lastPartDefForCS = aPartDef;
     437                if (PrefwdCS >0. &&  PreadjCS >0.) {
     438                        forward_CS_is_used = true;
     439                        LastCSCorrectionFactor = PrefwdCS/PreadjCS;
     440                }
     441                else {
     442                        forward_CS_is_used = false;
     443                        LastCSCorrectionFactor = 1.;
     444                       
     445                }
     446               
     447        }
     448        corr_fac =LastCSCorrectionFactor;
     449       
     450       
     451       
     452  } 
     453  else {
     454        forward_CS_is_used = false;
     455        LastCSCorrectionFactor = 1.;
     456  }
     457  fwd_TotCS=PrefwdCS;   
     458  fwd_is_used = forward_CS_is_used;
     459  return  corr_fac;
    290460}                                   
    291461///////////////////////////////////////////////////////
    … …  
    293463G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
    294464                                     const G4MaterialCutsCouple* aCouple, G4double step_length)
    295 { //G4double fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
    296  
    297   G4double corr_fac = 1.;
    298   if (consider_continuous_weight_correction) {
    299        
    300         G4double adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
    301         G4double PrefwdCS;
    302         PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
    303         G4double fwdCS = GetTotalForwardCS(aPartDef, (AfterStepEkin+PreStepEkin)/2.,aCouple);
    304         G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
    305         //if (aPartDef ==G4AdjointGamma::AdjointGamma()) G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
    306         /*if (adjCS >0 ) corr_fac = std::exp((PrefwdCS-fwdCS)*step_length);
    307         else corr_fac = std::exp(-fwdCS*step_length);*/
    308         corr_fac *=std::exp((adjCS-fwdCS)*step_length);
    309         corr_fac=std::max(corr_fac,1.e-6);
    310         corr_fac *=PreStepEkin/AfterStepEkin;
    311        
    312   }
    313   G4cout<<"Cont "<<corr_fac<<std::endl;
    314   G4cout<<"Ekin0 "<<PreStepEkin<<std::endl;
    315   G4cout<<"Ekin1 "<<AfterStepEkin<<std::endl;
    316   G4cout<<"step_length "<<step_length<<std::endl;
     465{  G4double corr_fac = 1.;
     466  //return corr_fac;
     467  //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
     468  G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
     469  G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
     470  if (!forward_CS_is_used || pre_adjCS ==0. ||  after_fwdCS==0.) {
     471        forward_CS_is_used=false;
     472        G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
     473        corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
     474        LastCSCorrectionFactor = 1.;
     475  }
     476  else {
     477        LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
     478  }     
     479       
     480
     481 
    317482  return corr_fac;
    318483}                                   
    319484///////////////////////////////////////////////////////
    320485//                                                           
    321 G4double G4AdjointCSManager::GetPostStepWeightCorrection(G4ParticleDefinition* , G4ParticleDefinition* ,
    322                                                           G4double ,G4double ,
    323                                                           const G4MaterialCutsCouple* )
    324 { G4double corr_fac = 1.;
    325   if (consider_poststep_weight_correction) {
    326         /*G4double fwdCS = GetTotalForwardCS(aSecondPartDef, EkinPrim,aCouple);
    327         G4double adjCS = GetTotalAdjointCS(aPrimPartDef, EkinPrim,aCouple);*/
    328         //G4double fwd1CS = GetTotalForwardCS(aPrimPartDef, EkinPrim,aCouple);
    329         //if (adjCS>0 && fwd1CS>0) adjCS = fwd1CS;
    330         //corr_fac =fwdCS*EkinSecond/adjCS/EkinPrim;
    331         //corr_fac = adjCS/fwdCS;
    332   }
    333   return corr_fac;
     486G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
     487{//return 1.;
     488 return  1./LastCSCorrectionFactor;
     489       
    334490}                                                         
    335491///////////////////////////////////////////////////////
    336492//
    337 double  G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
     493G4double  G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
    338494                                                         G4VEmAdjointModel* aModel,
    339495                                                         G4double PrimEnergy,
    340496                                                         G4double Tcut,
    341497                                                         G4bool IsScatProjToProjCase,
    342                                                          std::vector<double>& CS_Vs_Element)
     498                                                         std::vector<G4double>& CS_Vs_Element)
    343499{
    344500 
     501  G4double EminSec=0;
     502  G4double EmaxSec=0;
     503 
     504  if (IsScatProjToProjCase){
     505        EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
     506        EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
     507  }
     508  else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
     509        EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
     510        EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
     511  }
     512  if (EminSec >= EmaxSec) return 0.;
     513
     514
    345515  G4bool need_to_compute=false;
    346516  if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
    … …  
    381551        CS_Vs_Element.clear();
    382552        if (!aModel->GetUseMatrix()){
    383                 return aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase);
     553                CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
    384554                                         
    385555       
    … …  
    397567                                        CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
    398568                        G4double factor=0.;
    399                         for (size_t i=0;i<n_el;i++){
    400                                 size_t ind_el = aMaterial->GetElement(i)->GetIndex();
     569                        for (size_t i=0;i<n_el;i++){ //this could be computed only once
     570                                //size_t ind_el = aMaterial->GetElement(i)->GetIndex();
    401571                                factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
    402                                 G4AdjointCSMatrix* theCSMatrix;
    403                                 if (IsScatProjToProjCase){
    404                                         theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
    405                                 }
    406                                 else    theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
    407                                 //G4double CS =0.;
    408                                
    409                                 //G4cout<<CS<<std::endl;                       
    410                                
    411572                        }
    412573                        CS *=factor;
    … …  
    417578                        for (size_t i=0;i<n_el;i++){
    418579                                size_t ind_el = aMaterial->GetElement(i)->GetIndex();
    419                                 //G4cout<<aMaterial->GetName()<<std::endl;
     580                                //G4cout<<aMaterial->GetName()<<G4endl;
    420581                                G4AdjointCSMatrix* theCSMatrix;
    421582                                if (IsScatProjToProjCase){
    … …  
    426587                                if (PrimEnergy > Tlow)
    427588                                        CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
    428                                 //G4cout<<CS<<std::endl;                       
     589                                //G4cout<<CS<<G4endl;                   
    429590                                CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
    430591                        }
    … …  
    453614  G4double CS=0;
    454615  for (size_t i=0;i<CS_Vs_Element.size();i++){
    455         CS+=CS_Vs_Element[i];
    456   }
    457 
     616        CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
     617       
     618  }
    458619  return CS;
    459        
    460  
    461  
    462  
    463  
    464  
    465  
    466  
    467620}                                                       
    468621///////////////////////////////////////////////////////
    … …  
    473626                                                                   G4double Tcut,
    474627                                                                   G4bool IsScatProjToProjCase)
    475 { std::vector<double> CS_Vs_Element;
     628{ std::vector<G4double> CS_Vs_Element;
    476629  G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
    477630  G4double rand_var= G4UniformRand();
    … …  
    498651{
    499652 G4double TotalCS=0.;
    500 // G4ParticleDefinition* theDirPartDef = GetForwardParticleEquivalent(aPartDef);
     653
    501654 DefineCurrentMaterial(aCouple);
    502 /* size_t idx=-1;
    503  if (theDirPartDef->GetParticleName() == "gamma") idx = 0;
    504  else if (theDirPartDef->GetParticleName() == "e-") idx = 1;
    505  else if (theDirPartDef->GetParticleName() == "e+") idx = 2;
    506  
    507  //THe tCut computation is wrong this should be on Tcut per model the secondary determioming the Tcut
    508  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
    509  //G4cout<<aVec<<std::endl;
    510  G4double Tcut =(*aVec)[aCouple->GetIndex()];*/
    511  //G4cout<<"Tcut "<<Tcut<<std::endl;
    512  //G4cout<<(*aVec)[0]<<std::endl;
    513 // G4double Tcut =converters[idx]->Convert(Rcut,aCouple->GetMaterial());
    514  
    515  
    516  std::vector<double> CS_Vs_Element;
     655
     656 
     657 std::vector<G4double> CS_Vs_Element;
    517658 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
    518         /*G4ParticleDefinition* theDirSecondPartDef =
    519                         GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
    520659       
    521         */
    522        
    523        
    524660        G4double Tlow=0;
    525661        if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
    … …  
    527663                G4ParticleDefinition* theDirSecondPartDef =
    528664                        GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
    529                 G4int idx=-1;
     665                size_t idx=56;
    530666                if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
    531667                else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
    532668                else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
    533                 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
    534                 Tlow =(*aVec)[aCouple->GetIndex()];
     669                if (idx <56) {
     670                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
     671                        Tlow =(*aVec)[aCouple->GetIndex()];
     672                }       
    535673               
    536674       
    … …  
    539677        if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
    540678                if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
    541                         //G4cout<<"Yes1 before "<<std::endl;
    542679                        TotalCS += ComputeAdjointCS(currentMaterial,
    543680                                                       listOfAdjointEMModel[i],
    544                                                        Ekin, Tlow,true,CS_Vs_Element);
    545                         //G4cout<<"Yes1 "<<Ekin<<'\t'<<TotalCS<<std::endl;                             
     681                                                       Ekin, Tlow,true,CS_Vs_Element);                         
    546682                }
    547683                if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
    … …  
    549685                                                       listOfAdjointEMModel[i],
    550686                                                       Ekin, Tlow,false, CS_Vs_Element);
    551                                                        
    552                         //G4cout<<"Yes2 "<<TotalCS<<std::endl;
    553687                }
    554688               
    … …  
    563697std::vector<G4AdjointCSMatrix*>
    564698G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
    565                                                                         int nbin_pro_decade)
     699                                                                        G4int nbin_pro_decade)
    566700{
    567701  G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
    … …  
    577711   
    578712   
    579    
    580    
    581    
    582    
    583    
    584713   //Product to projectile backward scattering
    585714   //-----------------------------------------
    586715   G4double fE=std::pow(10.,1./nbin_pro_decade);
    587    G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
     716   G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    588717   G4double E1=EkinMin;
    589718   while (E1 <EkinMaxForProd){
    590719        E1=std::max(EkinMin,E2);
    591720        E1=std::min(EkinMaxForProd,E1);
    592         std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
     721        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
    593722        if (aMat.size()>=2) {
    594                 std::vector< G4double >* log_ESecVec=aMat[0];
    595                 std::vector< G4double >* log_CSVec=aMat[1];
     723                std::vector< double>* log_ESecVec=aMat[0];
     724                std::vector< double>* log_CSVec=aMat[1];
    596725                G4double log_adjointCS=log_CSVec->back();
    597726                //normalise CSVec such that it becomes a probability vector
    598                 /*for (size_t j=0;j<log_CSVec->size();j++) (*log_CSVec)[j]=(*log_CSVec)[j]-log_adjointCS;
    599                 (*log_CSVec)[0]=-90.;*/
    600        
    601        
    602                 for (size_t j=0;j<log_CSVec->size();j++) {
    603                         //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
    604                         if (j==0) (*log_CSVec)[j] = 0.;
    605                         else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
    606                         //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
     727                for (size_t j=0;j<log_CSVec->size();j++) {
     728                        if (j==0) (*log_CSVec)[j] = 0.;
     729                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
    607730                }       
    608                 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
     731                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
    609732                theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
    610733        }       
    … …  
    616739   //-----------------------------------------
    617740   
    618    E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
     741   E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    619742   E1=EkinMin;
    620743   while (E1 <EkinMaxForScat){
    621744        E1=std::max(EkinMin,E2);
    622745        E1=std::min(EkinMaxForScat,E1);
    623         std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
     746        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
    624747        if (aMat.size()>=2) {
    625                 std::vector< G4double >* log_ESecVec=aMat[0];
    626                 std::vector< G4double >* log_CSVec=aMat[1];
     748                std::vector< double>* log_ESecVec=aMat[0];
     749                std::vector< double>* log_CSVec=aMat[1];
    627750                G4double log_adjointCS=log_CSVec->back();
    628751                //normalise CSVec such that it becomes a probability vector
    629752                for (size_t j=0;j<log_CSVec->size();j++) {
    630                         //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
    631                         if (j==0) (*log_CSVec)[j] = 0.;
    632                         else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
    633                 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
     753                        if (j==0) (*log_CSVec)[j] = 0.;
     754                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
    634755                }       
    635                 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
     756                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
    636757                theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
    637758        }
    … …  
    641762   
    642763   
    643    
    644    
    645    
    646    
    647    
    648764  std::vector<G4AdjointCSMatrix*> res;
    649765  res.clear();
    650  
    651766  res.push_back(theCSMatForProdToProjBackwardScattering);
    652767  res.push_back(theCSMatForScatProjToProjBackwardScattering);
    653768 
    654769
    655 #ifdef TEST_MODE
     770/*
    656771  G4String file_name;
    657772  std::stringstream astream;
    … …  
    662777  theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
    663778 
    664   /*G4AdjointCSMatrix* aMat1 = new G4AdjointCSMatrix(false);
    665   G4AdjointCSMatrix* aMat2 = new G4AdjointCSMatrix(true);
    666  
    667   aMat1->Read(G4String("test_Z")+str_Z+"_1.txt");
    668   aMat2->Read(G4String("test_Z")+str_Z+"_2.txt");
    669   aMat1->Write(G4String("test_Z")+str_Z+"_11.txt");
    670   aMat2->Write(G4String("test_Z")+str_Z+"_22.txt"); */
    671 #endif
     779*/
     780
    672781 
    673782  return res;
    … …  
    702811   //-----------------------------------------
    703812   G4double fE=std::pow(10.,1./nbin_pro_decade);
    704    G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
     813   G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    705814   G4double E1=EkinMin;
    706815   while (E1 <EkinMaxForProd){
    707816        E1=std::max(EkinMin,E2);
    708817        E1=std::min(EkinMaxForProd,E1);
    709         std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
     818        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
    710819        if (aMat.size()>=2) {
    711                 std::vector< G4double >* log_ESecVec=aMat[0];
    712                 std::vector< G4double >* log_CSVec=aMat[1];
     820                std::vector< double>* log_ESecVec=aMat[0];
     821                std::vector< double>* log_CSVec=aMat[1];
    713822                G4double log_adjointCS=log_CSVec->back();
    714823       
    715824                //normalise CSVec such that it becomes a probability vector
    716825                for (size_t j=0;j<log_CSVec->size();j++) {
    717                         //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
     826                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
    718827                        if (j==0) (*log_CSVec)[j] = 0.;
    719828                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
    720                         //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
     829                        //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
    721830                }       
    722                 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
     831                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
    723832                theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
    724833        }       
    … …  
    733842   //-----------------------------------------
    734843   
    735    E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
     844   E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    736845   E1=EkinMin;
    737846   while (E1 <EkinMaxForScat){
    738847        E1=std::max(EkinMin,E2);
    739848        E1=std::min(EkinMaxForScat,E1);
    740         std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
     849        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
    741850        if (aMat.size()>=2) {
    742                 std::vector< G4double >* log_ESecVec=aMat[0];
    743                 std::vector< G4double >* log_CSVec=aMat[1];
     851                std::vector< double>* log_ESecVec=aMat[0];
     852                std::vector< double>* log_CSVec=aMat[1];
    744853                G4double log_adjointCS=log_CSVec->back();
    745854       
    746855                for (size_t j=0;j<log_CSVec->size();j++) {
    747                         //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
     856                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
    748857                        if (j==0) (*log_CSVec)[j] = 0.;
    749858                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
    750                 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
     859                //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
     860 
    751861                }       
    752                 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
     862                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
    753863       
    754864                theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
    … …  
    770880  res.push_back(theCSMatForScatProjToProjBackwardScattering);
    771881 
    772 #ifdef TEST_MODE
     882 /*
    773883  theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
    774884  theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
    775 #endif
     885*/
    776886
    777887
    … …  
    786896{
    787897 if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
    788  if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
     898 else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
     899 else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
     900 else if (theFwdPartDef ==theFwdIon) return theAdjIon;
     901 
    789902 return 0;     
    790903}
    … …  
    794907{
    795908 if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
    796  if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
     909 else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
     910 else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
     911 else if (theAdjPartDef == theAdjIon) return theFwdIon;
    797912 return 0;     
    798913}
    … …  
    805920    currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
    806921    currentMatIndex = couple->GetIndex();
    807     //G4cout<<"Index material "<<currentMatIndex<<std::endl;
     922    lastPartDefForCS =0;
     923    LastEkinForCS =0;
     924    LastCSCorrectionFactor =1.;
    808925  } 
    809926}
    810927
    811 
    812 
    813 ///////////////////////////////////////////////////////
    814 //
    815 double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
     928///////////////////////////////////////////////////////
     929//
     930void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
     931{
     932  if(aPartDef != currentParticleDef) {
     933 
     934        currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
     935        massRatio=1;
     936        if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
     937        currentParticleIndex=1000000;
     938        for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
     939                if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
     940        }       
     941         
     942  } 
     943}
     944
     945
     946
     947/////////////////////////////////////////////////////////////////////////////////////////////////
     948//
     949G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
    816950                                anAdjointCSMatrix,G4double Tcut)
    817951{
    818   std::vector< G4double > *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
     952  std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
    819953  if (theLogPrimEnergyVector->size() ==0){
    820         G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
    821         G4cout<<"The sampling procedure will be stopped."<<std::endl;
     954        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
     955        G4cout<<"The s"<<G4endl;
    822956        return 0.;
    823957       
    824958  }
    825   //G4cout<<"A prim/Tcut "<<aPrimEnergy<<'\t'<<Tcut<<std::endl;
    826959  G4double log_Tcut = std::log(Tcut);
    827960  G4double log_E =std::log(aPrimEnergy);
    … …  
    834967 
    835968  size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
    836   //G4cout<<"Prim energy "<<(*thePrimEnergyVector)[0]<<std::endl;
    837   //G4cout<<"Prim energy[ind]"<<(*thePrimEnergyVector)[ind]<<std::endl;
    838   //G4cout<<"Prim energy ind"<<ind<<std::endl;
    839  
    840969  G4double aLogPrimEnergy1,aLogPrimEnergy2;
    841970  G4double aLogCS1,aLogCS2;
    842971  G4double log01,log02;
    843   std::vector< G4double>* aLogSecondEnergyVector1 =0;
    844   std::vector< G4double>* aLogSecondEnergyVector2  =0;
    845   std::vector< G4double>* aLogProbVector1=0;
    846   std::vector< G4double>* aLogProbVector2=0;
     972  std::vector< double>* aLogSecondEnergyVector1 =0;
     973  std::vector< double>* aLogSecondEnergyVector2  =0;
     974  std::vector< double>* aLogProbVector1=0;
     975  std::vector< double>* aLogProbVector2=0;
    847976  std::vector< size_t>* aLogProbVectorIndex1=0;
    848977  std::vector< size_t>* aLogProbVectorIndex2=0;
    … …  
    851980  anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
    852981  anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
    853   //G4cout<<"aSecondEnergyVector1.size() "<<aSecondEnergyVector1->size()<<std::endl;
    854   //G4cout<<aSecondEnergyVector1<<std::endl;
    855   //G4cout<<"aSecondEnergyVector2.size() "<<aSecondEnergyVector2->size()<<std::endl;
    856982  if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
    857983        G4double log_minimum_prob1, log_minimum_prob2;
    858        
    859         //G4cout<<aSecondEnergyVector1->size()<<std::endl;
    860984        log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
    861985        log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
    862         //G4cout<<"minimum_prob1 "<< std::exp(log_minimum_prob1)<<std::endl;
    863         //G4cout<<"minimum_prob2 "<< std::exp(log_minimum_prob2)<<std::endl;
    864         //G4cout<<"Tcut "<<std::endl;
    865986        aLogCS1+= log_minimum_prob1;
    866987        aLogCS2+= log_minimum_prob2;

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointCSMatrix.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
     29
    2630#include "G4AdjointCSMatrix.hh"
    2731#include <iomanip>
    2832#include <fstream>
    29 
    3033#include "G4AdjointInterpolator.hh"
    31 
    3234///////////////////////////////////////////////////////
    3335//
    … …  
    6466///////////////////////////////////////////////////////
    6567//
    66  void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< G4double>* aLogSecondEnergyVector,
    67                                                                        std::vector< G4double>* aLogProbVector,size_t n_pro_decade){
     68 void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< double>* aLogSecondEnergyVector,
     69                                                                       std::vector< double>* aLogProbVector,size_t n_pro_decade){
    6870       
    6971        G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
    70         //Add this time we consider that the energy are given monotically
    7172       
     73        //At this time we consider that the energy is increasing monotically
    7274        theLogPrimEnergyVector.push_back(aLogPrimEnergy);
    7375        theLogCrossSectionVector.push_back(aLogCS);
    7476        theLogSecondEnergyMatrix.push_back(aLogSecondEnergyVector);
    75         //G4cout<<"Test Add Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
    76         //G4cout<<theSecondEnergyMatrix.size()<<std::endl;
    7777        theLogProbMatrix.push_back(aLogProbVector);
    78         //G4cout<<"Test Add Data 1 "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
    79         //G4cout<<theSecondEnergyMatrix.size()<<std::endl;
     78       
    8079        std::vector< size_t>* aLogProbVectorIndex = 0;
    8180        dlog =0;
     81       
    8282        if (n_pro_decade > 0 && aLogProbVector->size()>0) {
    8383                aLogProbVectorIndex = new std::vector< size_t>();
    … …  
    8585                G4double log_val = int(std::min((*aLogProbVector)[0],aLogProbVector->back())/dlog)*dlog;
    8686                log0Vector.push_back(log_val);
     87               
    8788                while(log_val<0.) {
    8889                        aLogProbVectorIndex->push_back(theInterpolator->FindPosition(log_val,(*aLogProbVector)));
    … …  
    102103///////////////////////////////////////////////////////
    103104//
    104 bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< G4double>*& aLogSecondEnergyVector,
    105                                                                        std::vector< G4double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
     105G4bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< double>*& aLogSecondEnergyVector,
     106                                                                       std::vector< double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
    106107{       if (i>= nb_of_PrimEnergy) return false;
    107         //G4cout<<"Test Get Data "<<std::endl;
     108        //G4cout<<"Test Get Data "<<G4endl;
    108109        aLogPrimEnergy = theLogPrimEnergyVector[i];
    109110        aLogCS = theLogCrossSectionVector[i];
    110111        aLogSecondEnergyVector = theLogSecondEnergyMatrix[i];
    111         //G4cout<<"Test Get Data "<<this<<'\t'<<theSecondEnergyMatrix[i]->size()<<std::endl;
    112         //G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
    113         //G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector<<std::endl;
    114112        aLogProbVector = theLogProbMatrix[i];
    115113        aLogProbVectorIndex = theLogProbMatrixIndex[i];
    116114        log0=log0Vector[i];
    117         //G4cout<<"Test Get Data 1 "<<this<<'\t'<<theProbMatrix[i]->size()<<std::endl;
    118         //G4cout<<"Test Get Data 1 "<<this<<'\t'<<aProbVector->size()<<std::endl;
    119         //G4cout<<"Test Get Data 1 "<<this<<'\t'<<aLogProbVectorIndex<<std::endl;
    120115        return true;
    121116       
    … …  
    127122        FileOutput<<std::setiosflags(std::ios::scientific);
    128123        FileOutput<<std::setprecision(6);
    129         FileOutput<<theLogPrimEnergyVector.size()<<std::endl;
     124        FileOutput<<theLogPrimEnergyVector.size()<<G4endl;
    130125        for (size_t i=0;i<theLogPrimEnergyVector.size();i++){
    131                 FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<std::endl;
     126                FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<G4endl;
    132127                size_t j1=0;
    133                 FileOutput<<theLogSecondEnergyMatrix[i]->size()<<std::endl;
     128                FileOutput<<theLogSecondEnergyMatrix[i]->size()<<G4endl;
    134129                for (size_t j=0;j<theLogSecondEnergyMatrix[i]->size();j++){
    135130                        FileOutput<<std::exp((*theLogSecondEnergyMatrix[i])[j]);
    … …  
    137132                        if (j1<10) FileOutput<<'\t';
    138133                        else  {
    139                                 FileOutput<<std::endl;
     134                                FileOutput<<G4endl;
    140135                                j1=0;
    141136                        }       
    142137                }
    143                 if (j1>0) FileOutput<<std::endl;
     138                if (j1>0) FileOutput<<G4endl;
    144139                j1=0;
    145                 FileOutput<<theLogProbMatrix[i]->size()<<std::endl;
     140                FileOutput<<theLogProbMatrix[i]->size()<<G4endl;
    146141                for (size_t j=0;j<theLogProbMatrix[i]->size();j++){
    147142                        FileOutput<<std::exp((*theLogProbMatrix[i])[j]);
    … …  
    149144                        if (j1<10) FileOutput<<'\t';
    150145                        else  {
    151                                 FileOutput<<std::endl;
     146                                FileOutput<<G4endl;
    152147                                j1=0;
    153148                        }       
    154149                }
    155                 if (j1>0) FileOutput<<std::endl;
     150                if (j1>0) FileOutput<<G4endl;
    156151               
    157152               
    … …  
    177172                theLogCrossSectionVector.push_back(CS);
    178173                FileOutput>>n2;
    179                 theLogSecondEnergyMatrix.push_back(new std::vector<double>());
    180                 theLogProbMatrix.push_back(new std::vector<double>());
     174                theLogSecondEnergyMatrix.push_back(new std::vector<G4double>());
     175                theLogProbMatrix.push_back(new std::vector<G4double>());
    181176               
    182177                for (size_t j=0; j<n2;j++){

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointComptonModel.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointComptonModel.hh"
    2730#include "G4AdjointCSManager.hh"
    … …  
    3437#include "G4AdjointGamma.hh"
    3538#include "G4Gamma.hh"
     39#include "G4KleinNishinaCompton.hh"
    3640
    3741
    … …  
    4448  SetUseMatrix(true);
    4549  SetUseMatrixPerElement(true);
    46   SetIsIonisation(false);
    4750  SetUseOnlyOneMatrixForAllElements(true);
    4851  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
    … …  
    5053  theDirectPrimaryPartDef=G4Gamma::Gamma();
    5154  second_part_of_same_type=false;
     55  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
    5256 
    5357}
    … …  
    6266                       G4ParticleChange* fParticleChange)
    6367{
    64    
     68   if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
    6569   
    6670   //A recall of the compton scattering law is
    … …  
    7175   
    7276  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
    73   //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
    74   size_t ind= 0;
    75  
    76  
    77  
    78  
    79  //Elastic inverse scattering //not correct in all the cases
    80  //---------------------------------------------------------
    81  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
    82  
    83  //G4cout<<adjointPrimKinEnergy<<std::endl;
    84  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
     77  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
     78  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
    8579        return;
    86  }
     80  }
     81 
    8782 
    8883 //Sample secondary energy
    8984 //-----------------------
    9085 G4double gammaE1;
    91 
    92  gammaE1 = SampleAdjSecEnergyFromCSMatrix(ind,
    93                                                   adjointPrimKinEnergy,
     86 gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
    9487                                                  IsScatProjToProjCase);
    9588 
    … …  
    108101 //Cos th
    109102 //-------
    110 // G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<std::endl;
     103// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl;
    111104 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
    112105 if (!IsScatProjToProjCase) {
    … …  
    116109 G4double sin_th = 0.;
    117110 if (std::abs(cos_th)>1){
    118         //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<std::endl;
     111        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
    119112        if (cos_th>0) {
    120113                cos_th=1.;
    … …  
    136129 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
    137130 gammaMomentum1.rotateUz(dir_parallel);
    138 // G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<std::endl;
     131// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl;
    139132 
    140133 
    141134 //It is important to correct the weight of particles before adding the secondary
    142135 //------------------------------------------------------------------------------
    143  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,gammaE1);
    144  
    145  if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
     136 CorrectPostStepWeight(fParticleChange,
     137                        aTrack.GetWeight(),
     138                        adjointPrimKinEnergy,
     139                        gammaE1,
     140                        IsScatProjToProjCase);
     141 
     142 if (!IsScatProjToProjCase){ //kill the primary and add a secondary
    146143        fParticleChange->ProposeTrackStatus(fStopAndKill);
    147144        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
    148         //G4cout<<"gamma0Momentum "<<gamma0Momentum<<std::endl;
     145        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
    149146 }
    150147 else {
    … …  
    154151 
    155152       
    156 }                       
     153}
     154////////////////////////////////////////////////////////////////////////////////
     155//
     156void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
     157                       G4bool IsScatProjToProjCase,
     158                       G4ParticleChange* fParticleChange)
     159{
     160
     161 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
     162 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
     163 
     164 
     165 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
     166
     167 
     168 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
     169        return;
     170 }
     171 
     172 
     173 G4double diffCSUsed=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
     174 G4double gammaE1=0.;
     175 G4double gammaE2=0.;
     176 if (!IsScatProjToProjCase){
     177       
     178        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
     179        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
     180        if (Emin>=Emax) return;
     181        G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
     182        G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
     183        gammaE1=adjointPrimKinEnergy/(1.-f1*pow(f2,G4UniformRand()));;
     184        gammaE2=gammaE1-adjointPrimKinEnergy;
     185        diffCSUsed= diffCSUsed*(1.+2.*log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
     186       
     187       
     188 }
     189 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
     190        G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
     191        if (Emin>=Emax) return;
     192        gammaE2 =adjointPrimKinEnergy;
     193        gammaE1=Emin*pow(Emax/Emin,G4UniformRand());
     194        diffCSUsed= diffCSUsed/gammaE1;
     195       
     196 }
     197 
     198 
     199 
     200                                                       
     201 //Weight correction
     202 //-----------------------
     203 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
     204 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
     205
     206 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
     207 //one consistent with the direct model
     208 
     209 
     210 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
     211 if (diffCS >0)  diffCS /=G4direct_CS;  // here we have the normalised diffCS
     212 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
     213 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
     214 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;                                 
     215 
     216 w_corr*=diffCS/diffCSUsed;
     217           
     218 G4double new_weight = aTrack.GetWeight()*w_corr;
     219 fParticleChange->SetParentWeightByProcess(false);
     220 fParticleChange->SetSecondaryWeightByProcess(false);
     221 fParticleChange->ProposeParentWeight(new_weight);
     222 
     223 
     224 
     225 //Cos th
     226 //-------
     227
     228 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
     229 if (!IsScatProjToProjCase) {
     230        G4double p_elec=theAdjointPrimary->GetTotalMomentum();
     231        cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
     232 }
     233 G4double sin_th = 0.;
     234 if (std::abs(cos_th)>1){
     235        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
     236        if (cos_th>0) {
     237                cos_th=1.;
     238        }
     239        else    cos_th=-1.;
     240        sin_th=0.;
     241 }
     242 else  sin_th = std::sqrt(1.-cos_th*cos_th);
     243
     244 
     245 
     246 
     247 //gamma0 momentum
     248 //--------------------
     249
     250 
     251 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
     252 G4double phi =G4UniformRand()*2.*3.1415926;
     253 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
     254 gammaMomentum1.rotateUz(dir_parallel);
     255 
     256 
     257
     258 
     259 if (!IsScatProjToProjCase){ //kill the primary and add a secondary
     260        fParticleChange->ProposeTrackStatus(fStopAndKill);
     261        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
     262        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
     263 }
     264 else {
     265        fParticleChange->ProposeEnergy(gammaE1);
     266        fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
     267 }
     268 
     269 
     270 
     271}
     272
     273                       
    157274////////////////////////////////////////////////////////////////////////////////
    158275//
    … …  
    179296                                      G4double )
    180297{ //Based on Klein Nishina formula
    181  //* In the forward case (see G4KleinNishinaModel)  the  cross section is parametrised while the secondaries are sampled from the
     298 // In the forward case (see G4KleinNishinaModel)  the  cross section is parametrised while
     299 // the secondaries are sampled from the
    182300 // Klein Nishida differential cross section
    183  // The used diffrential cross section here is therefore the cross section multiplied by the normalidsed differential Klein Nishida cross section
     301 // The used diffrential cross section here is therefore the cross section multiplied by the normalised
     302 //differential Klein Nishida cross section
    184303 
    185304 
    … …  
    191310 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
    192311 if (gamEnergy1 >gamEnergy1_max ||  gamEnergy1<gamEnergy1_min) {
    193         /*G4cout<<"the differential CS is null"<<std::endl;
    194         G4cout<<gamEnergy0<<std::endl;
    195         G4cout<<gamEnergy1<<std::endl;
    196         G4cout<<gamEnergy1_min<<std::endl;*/
     312        /*G4cout<<"the differential CS is null"<<G4endl;
     313        G4cout<<gamEnergy0<<G4endl;
     314        G4cout<<gamEnergy1<<G4endl;
     315        G4cout<<gamEnergy1_min<<G4endl;*/
    197316        return 0.;
    198317 }
    … …  
    222341 //-------------------------------
    223342 
    224  G4double G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
     343 G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
    225344                                             gamEnergy0,
    226345                                             Z, 0., 0.,0.);
    227346 
    228347 dCS_dE1 *= G4direct_CS/CS;
    229 /* G4cout<<"the differential CS is not null"<<std::endl;
    230  G4cout<<gamEnergy0<<std::endl;
    231  G4cout<<gamEnergy1<<std::endl;*/
     348/* G4cout<<"the differential CS is not null"<<G4endl;
     349 G4cout<<gamEnergy0<<G4endl;
     350 G4cout<<gamEnergy1<<G4endl;*/
    232351 
    233352 return dCS_dE1;
    … …  
    235354
    236355}
     356                                     
    237357////////////////////////////////////////////////////////////////////////////////
    238358//
    … …  
    251371  return  emin;
    252372}
     373////////////////////////////////////////////////////////////////////////////////
     374//
     375G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
     376                                             G4double primEnergy,
     377                                             G4bool IsScatProjToProjCase)
     378{
     379  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
     380  DefineCurrentMaterial(aCouple);
     381 
     382 
     383  G4double Cross=0.;
     384  G4double Emax_proj =0.;
     385  G4double Emin_proj =0.;
     386  if (!IsScatProjToProjCase ){
     387        Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
     388        Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
     389        if (Emax_proj>Emin_proj ){
     390                 Cross= log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
     391                                                                *(1.+2.*log(1.+electron_mass_c2/primEnergy));
     392        }       
     393  }
     394  else {
     395        Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
     396        Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
     397        if (Emax_proj>Emin_proj) {
     398                Cross = log(Emax_proj/Emin_proj);
     399                //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
     400        }
     401       
     402       
     403  }
     404 
     405  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
     406  lastCS=Cross;
     407  return Cross;
     408}

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointInterpolator.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointCSMatrix.hh"
    2730#include "G4AdjointInterpolator.hh"
    … …  
    6164G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
    6265{ G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
    63   //G4cout<<"Linear "<<res<<std::endl;
     66  //G4cout<<"Linear "<<res<<G4endl;
    6467  return res;   
    6568}
    … …  
    6972{ if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
    7073  G4double B=std::log(y2/y1)/std::log(x2/x1);
    71   //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<std::endl;
     74  //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<G4endl;
    7275  G4double A=y1/std::pow(x1,B);
    7376  G4double res=A*std::pow(x,B);
    74  // G4cout<<"Log "<<res<<std::endl;
     77 // G4cout<<"Log "<<res<<G4endl;
    7578  return res;
    7679}
    … …  
    98101  }
    99102  else {
    100         //G4cout<<"The interpolation method that you invoked does not exist!"<<std::endl;
     103        //G4cout<<"The interpolation method that you invoked does not exist!"<<G4endl;
    101104        return -1111111111.;
    102105  }
    … …  
    104107///////////////////////////////////////////////////////
    105108//
    106 size_t  G4AdjointInterpolator::FindPosition(G4double& x,std::vector<double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
     109size_t  G4AdjointInterpolator::FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
    107110{  //most rapid nethod could be used probably
    108    //It is important to put std::vector<double>& such that the vector itself is used and not a copy
     111   //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
    109112 
    110113 
    … …  
    149152///////////////////////////////////////////////////////
    150153//
    151 size_t  G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<double>& log_x_vec) //only valid if x_vec is monotically increasing
     154size_t  G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec) //only valid if x_vec is monotically increasing
    152155{  //most rapid nethod could be used probably
    153    //It is important to put std::vector<double>& such that the vector itself is used and not a copy
    154  
     156   //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
     157  return FindPosition(log_x, log_x_vec);
    155158  if (log_x_vec.size()>3){
    156159        size_t ind=0;
    … …  
    170173///////////////////////////////////////////////////////
    171174//
    172 G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,G4String InterPolMethod)
     175G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
    173176{ size_t i=FindPosition(x,x_vec);
    174   //G4cout<<i<<std::endl;
    175   //G4cout<<x<<std::endl;
    176   //G4cout<<x_vec[i]<<std::endl;
     177  //G4cout<<i<<G4endl;
     178  //G4cout<<x<<G4endl;
     179  //G4cout<<x_vec[i]<<G4endl;
    177180  return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
    178181}
    … …  
    180183///////////////////////////////////////////////////////
    181184//
    182 G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,
     185G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
    183186                                            std::vector<size_t>& index_vec,G4double x0, G4double dx) //only linear interpolation possible
    184187{ size_t ind=0;
    … …  
    202205///////////////////////////////////////////////////////
    203206//
    204 G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<double>& log_x_vec,std::vector<double>& log_y_vec)
     207G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
    205208{ //size_t i=0;
    206209  size_t i=FindPositionForLogVector(log_x,log_x_vec);
    207   /*G4cout<<"In interpolate "<<std::endl;
    208   G4cout<<i<<std::endl;
    209   G4cout<<log_x<<std::endl;
    210   G4cout<<log_x_vec[i]<<std::endl;
    211   G4cout<<log_x_vec[i+1]<<std::endl;
    212   G4cout<<log_y_vec[i]<<std::endl;
    213   G4cout<<log_y_vec[i+1]<<std::endl;*/
     210  /*G4cout<<"In interpolate "<<G4endl;
     211  G4cout<<i<<G4endl;
     212  G4cout<<log_x<<G4endl;
     213  G4cout<<log_x_vec[i]<<G4endl;
     214  G4cout<<log_x_vec[i+1]<<G4endl;
     215  G4cout<<log_y_vec[i]<<G4endl;
     216  G4cout<<log_y_vec[i+1]<<G4endl;*/
    214217 
    215218  G4double log_y=LinearInterpolation(log_x,log_x_vec[i],log_x_vec[i+1],log_y_vec[i],log_y_vec[i+1]);

  • trunk/source/processes/electromagnetic/adjoint/src/G4AdjointPhotoElectricModel.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4AdjointPhotoElectricModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4AdjointPhotoElectricModel.hh"
    2730#include "G4AdjointCSManager.hh"
    … …  
    3538#include "G4AdjointGamma.hh"
    3639
     40
    3741////////////////////////////////////////////////////////////////////////////////
    3842//
    … …  
    4145
    4246{ SetUseMatrix(false);
     47  SetApplyCutInRange(false);
    4348  current_eEnergy =0.;
    4449  totAdjointCS=0.;
     50  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
     51  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
     52  theDirectPrimaryPartDef=G4Gamma::Gamma();
     53  second_part_of_same_type=false;
     54  theDirectPEEffectModel = new G4PEEffectModel();
     55 
    4556}
    4657////////////////////////////////////////////////////////////////////////////////
    … …  
    5768
    5869  //Compute the totAdjointCS vectors if not already done for the current couple and electron energy
     70  //-----------------------------------------------------------------------------------------------
    5971  const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
    6072  const G4DynamicParticle* aDynPart =  aTrack.GetDynamicParticle() ;
    6173  G4double electronEnergy = aDynPart->GetKineticEnergy();
    6274  G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
    63   totAdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
     75  pre_step_AdjointCS = totAdjointCS; //The last computed CS was  at pre step point
     76  G4double adjCS;
     77  adjCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
     78  post_step_AdjointCS = totAdjointCS;
    6479                               
    65  
    66   //Sample gamma energy
    67   //-------------
    68         /////////////////////////////////////////////////////////////////////////////////
    69 //      Module:         G4ContinuousGainOfEnergy.hh
    70 //      Author:         L. Desorgher
    71 //      Date:           1 September 2007
    72 //      Organisation:   SpaceIT GmbH
    73 //      Customer:       ESA/ESTEC
    74 /////////////////////////////////////////////////////////////////////////////////
    75 //
    76 // CHANGE HISTORY
    77 // --------------
    78 //      ChangeHistory:
    79 //              1 September 2007 creation by L. Desorgher               
    80 //
    81 //-------------------------------------------------------------
    82 //      Documentation:
    83 //              Modell for the adjoint compton scattering
    84 //
     80
     81
    8582
    8683  //Sample element
    8784  //-------------
    8885   const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
    89    const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
    9086   size_t nelm =  currentMaterial->GetNumberOfElements();
    91    G4double rand_CS= totAdjointCS*G4UniformRand();
     87   G4double rand_CS= G4UniformRand()*xsec[nelm-1];
    9288   for (index_element=0; index_element<nelm-1; index_element++){
    9389        if (rand_CS<xsec[index_element]) break;
    … …  
    9692   //Sample shell and binding energy
    9793   //-------------
    98    rand_CS= totAdjointCS*G4UniformRand()/theAtomNumDensityVector[index_element];
    9994   G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
     95   rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand();
    10096   G4int i  = 0; 
    10197   for (i=0; i<nShells-1; i++){
    … …  
    113109   G4double gamma   = 1. + electronEnergy/electron_mass_c2;
    114110   if (gamma <= 5.) {
    115         G4double beta  = std::sqrt(gamma*gamma-1.)/gamma;
     111        G4double beta  = sqrt(gamma*gamma-1.)/gamma;
    116112        G4double b     = 0.5*gamma*(gamma-1.)*(gamma-2);
    117113
    … …  
    131127 
    132128     
    133   G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
     129  G4double sin_theta = sqrt(1.-cos_theta*cos_theta);
    134130  G4double Phi     = twopi * G4UniformRand();
    135   G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
     131  G4double dirx = sin_theta*cos(Phi),diry = sin_theta*sin(Phi),dirz = cos_theta;
    136132  G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz);
    137133  adjoint_gammaDirection.rotateUz(electronDirection);
    … …  
    141137  //Weight correction
    142138 //-----------------------                                         
    143   CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy);       
     139  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase); 
    144140 
    145141 
    … …  
    149145  G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
    150146                       G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
     147 
     148 
     149 
     150 
     151 
    151152  fParticleChange->ProposeTrackStatus(fStopAndKill);
    152   fParticleChange->AddSecondary(anAdjointGamma);
    153      
     153  fParticleChange->AddSecondary(anAdjointGamma);   
    154154       
    155155
    156156 
    157157
     158}
     159
     160////////////////////////////////////////////////////////////////////////////////
     161//
     162void G4AdjointPhotoElectricModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
     163                                                            G4double old_weight, 
     164                                                            G4double adjointPrimKinEnergy,
     165                                                            G4double projectileKinEnergy ,
     166                                                            G4bool  )
     167{
     168 G4double new_weight=old_weight;
     169
     170 G4double w_corr =G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection()/factorCSBiasing;
     171 w_corr*=post_step_AdjointCS/pre_step_AdjointCS;
     172 new_weight*=w_corr;
     173 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
     174 fParticleChange->SetParentWeightByProcess(false);
     175 fParticleChange->SetSecondaryWeightByProcess(false);
     176 fParticleChange->ProposeParentWeight(new_weight);     
    158177}       
    159178
    … …  
    164183                                G4double electronEnergy,
    165184                                G4bool IsScatProjToProjCase)
    166 { if (IsScatProjToProjCase) return 0.;
     185{
     186 
     187
     188  if (IsScatProjToProjCase)  return 0.;
     189
     190       
    167191  if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
    168192        totAdjointCS = 0.;
    169193        DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
    170194        const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
    171         const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
     195        const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
    172196        size_t nelm =  currentMaterial->GetNumberOfElements();
    173197        for (index_element=0;index_element<nelm;index_element++){
    … …  
    176200                xsec[index_element] = totAdjointCS;
    177201        }
    178   }     
    179   return totAdjointCS;
    180  
     202
     203        totBiasedAdjointCS=std::min(totAdjointCS,0.01);
     204//      totBiasedAdjointCS=totAdjointCS;
     205        factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
     206        lastCS=totBiasedAdjointCS;
     207       
     208       
     209  }
     210  return totBiasedAdjointCS;
     211
    181212 
    182213}
    … …  
    188219  G4int nShells = anElement->GetNbOfAtomicShells();
    189220  G4double Z= anElement->GetZ();
    190   G4double N= anElement->GetN();
    191221  G4int i  = 0; 
    192222  G4double B0=anElement->GetAtomicShell(0);
    193223  G4double gammaEnergy = electronEnergy+B0;
    194   G4double adjointCS = theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,N,0.,0.)*electronEnergy/gammaEnergy;
     224  G4double CS= theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
     225  G4double adjointCS =0.;
     226  if (CS >0) adjointCS += CS/gammaEnergy;
    195227  shell_prob[index_element][0] = adjointCS;                                         
    196228  for (i=1;i<nShells;i++){
    197         //G4cout<<i<<std::endl;
     229        //G4cout<<i<<G4endl;
    198230        G4double Bi_= anElement->GetAtomicShell(i-1);
    199231        G4double Bi = anElement->GetAtomicShell(i);
    200         //G4cout<<Bi_<<'\t'<<Bi<<std::endl;
     232        //G4cout<<Bi_<<'\t'<<Bi<<G4endl;
    201233        if (electronEnergy <Bi_-Bi) {
    202234                gammaEnergy = electronEnergy+Bi;
    203                 adjointCS +=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,anElement->GetZ(),N,0.,0.)*electronEnergy/gammaEnergy;
     235               
     236                CS=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
     237                if (CS>0) adjointCS +=CS/gammaEnergy;
    204238        }
    205239        shell_prob[index_element][i] = adjointCS;       
    206240 
    207241  }
    208  
     242  adjointCS*=electronEnergy;
    209243  return adjointCS;
    210244 

  • trunk/source/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4ContinuousGainOfEnergy.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4ContinuousGainOfEnergy.hh"
    2730#include "G4Step.hh"
    … …  
    3134#include "G4VParticleChange.hh"
    3235#include "G4UnitsTable.hh"
     36#include "G4AdjointCSManager.hh"
     37#include "G4LossTableManager.hh"
    3338
    3439
    … …  
    3843  G4ProcessType type): G4VContinuousProcess(name, type)
    3944{
     45
    4046
    4147  linLossLimit=0.05;
    … …  
    4450  is_integral = false;
    4551 
     52  //Will be properly set in SetDirectParticle()
     53  IsIon=false;
     54  massRatio =1.;
     55  chargeSqRatio=1.;
     56  preStepChargeSqRatio=1.;
     57
     58 
     59
     60 
     61 
    4662}
    4763
    … …  
    7086}
    7187
    72 
    73 
     88///////////////////////////////////////////////////////
     89//
     90void  G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
     91{theDirectPartDef=p;
     92 if (theDirectPartDef->GetParticleType()== "nucleus") {
     93         IsIon=true;
     94         massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass();
     95         G4double q=theDirectPartDef->GetPDGCharge();
     96         chargeSqRatio=q*q;
     97       
     98         
     99 }
     100 
     101}
    74102
    75103///////////////////////////////////////////////////////
    … …  
    80108{
    81109   
     110  //Caution in this method the  step length should be the true step length
     111  // A problem is that this is compute by the multiple scattering that does not know the energy at the end of the adjoint step. This energy is used during the
     112  //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
     113  //
     114 
     115 
     116 
    82117  aParticleChange.Initialize(track);
    83118 
    … …  
    86121  G4double length = step.GetStepLength();
    87122  G4double degain  = 0.0;
    88  
    89  
    90  
    91 
     123 
     124 
     125 
    92126  // Compute this for weight change after continuous energy loss
    93127  //-------------------------------------------------------------
    94   G4double DEDX_before =
    95                 theDirectEnergyLossProcess
    96                                 ->GetDEDX(preStepKinEnergy, currentCouple);
    97  
    98  
     128  G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple);
     129   
     130 
    99131 
    100132  // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
    … …  
    103135  G4DynamicParticle* dynParticle = new G4DynamicParticle();
    104136  *dynParticle = *(track.GetDynamicParticle());
    105   G4double Tkin = dynParticle->GetKineticEnergy();
    106  
     137  dynParticle->SetDefinition(theDirectPartDef);
     138  G4double Tkin = dynParticle->GetKineticEnergy();
     139  G4double Tkin1=Tkin*0.001;
     140
    107141  size_t n=1;
    108142  if (is_integral ) n=10;
     143  n=1;
    109144  G4double dlength= length/n;
    110145  for (size_t i=0;i<n;i++) {
     146        G4double factor_dE=1.;
     147        if (Tkin != preStepKinEnergy && IsIon) {
     148                chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
     149                theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
     150       
     151        }
     152 
    111153        G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple);
    112154        if( dlength <= linLossLimit * r ) {
    113155                degain = DEDX_before*dlength;
     156                G4double degain1 = dlength*theDirectEnergyLossProcess->GetDEDX(Tkin1, currentCouple);
     157                factor_dE=1.+(degain1-degain)/(Tkin1-Tkin);
    114158        }
    115159        else {
    116                 G4double x = r + length;
    117                 degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
     160                G4double x = r + dlength;
     161                //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
     162                G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
     163                if (IsIon){
     164                        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
     165                        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
     166                        G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
     167                        while (std::abs(x-x1)>0.01*x) {
     168                                E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
     169                                chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
     170                                theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
     171                                x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
     172                       
     173                        }
     174                }
     175                G4double r1 = theDirectEnergyLossProcess->GetRange(Tkin1, currentCouple);
     176                G4double x1 = r1 + dlength;
     177                G4double E1 = theDirectEnergyLossProcess->GetKineticEnergy(x1,currentCouple);
     178                factor_dE=(E1-E)/(Tkin1-Tkin);
     179                degain=E-Tkin; 
     180               
     181               
     182               
    118183        }
    119         G4VEmModel* currentModel =  theDirectEnergyLossProcess->SelectModelForMaterial(Tkin+degain,currentMaterialIndex);
     184        //G4cout<<degain<<G4endl;
    120185        G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
    121186        tmax = std::min(tmax,currentTcut);
    122  
    123  
     187       
     188       
     189        dynParticle->SetKineticEnergy(Tkin+degain);
     190
     191        // Corrections, which cannot be tabulated for ions
     192        //----------------------------------------
     193        G4double esecdep=0;//not used in most models
     194        currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
    124195
    125196        // Sample fluctuations
    126197        //-------------------
     198       
    127199       
    128200        G4double deltaE =0.;
    129201        if (lossFluctuationFlag ) {
    130202                deltaE = currentModel->GetModelOfFluctuations()->
    131                                                 SampleFluctuations(currentMaterial,dynParticle,tmax,length,degain)-degain;
     203                                                SampleFluctuations(currentMaterial,dynParticle,tmax,dlength,degain)-degain;
    132204        }
    133         Tkin+=degain+deltaE;
     205       
     206        G4double egain=degain+deltaE;
     207        if (egain <=0) egain=degain;
     208        Tkin+=egain;
    134209        dynParticle->SetKineticEnergy(Tkin);
    135        
     210 }
     211 
     212 
     213 
     214
     215 
     216  delete dynParticle;
     217 
     218  if (IsIon){
     219        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
     220        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
     221               
    136222  }
    137  
    138  
    139 
    140   // Corrections, which cannot be tabulated
    141   // probably  this should be also changed
    142   // at this time it does nothing so we can leave it
    143   //CorrectionsAlongStep(currentCouple, dynParticle, egain, length);
    144  
    145   delete dynParticle;
    146  
    147223 
    148224  G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
    149   G4double weight_correction=DEDX_after/DEDX_before; //probably not needed
    150   weight_correction=1.;
     225 
     226 
     227  G4double weight_correction=DEDX_after/DEDX_before;
     228 
    151229 
    152230  aParticleChange.ProposeEnergy(Tkin);
    … …  
    154232  //we still need to register in the particleChange the modification of the weight of the particle
    155233  G4double new_weight=weight_correction*track.GetWeight();
    156   aParticleChange.SetParentWeightByProcess(true);
     234  aParticleChange.SetParentWeightByProcess(false);
    157235  aParticleChange.ProposeParentWeight(new_weight);
    158236 
    … …  
    168246  lossFluctuationFlag = val;
    169247}
     248///////////////////////////////////////////////////////
     249//
     250
     251
     252
     253G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
     254                G4double , G4double , G4double& )
     255{
     256  G4double x = DBL_MAX;
     257  x=.1*mm;
     258 
     259 
     260  DefineMaterial(track.GetMaterialCutsCouple());
     261 
     262  preStepKinEnergy = track.GetKineticEnergy();
     263  preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
     264  currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
     265  G4double emax_model=currentModel->HighEnergyLimit();
     266  if (IsIon) {
     267        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
     268        preStepChargeSqRatio = chargeSqRatio;
     269        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
     270  }
     271 
     272 
     273  G4double maxE =1.1*preStepKinEnergy;
     274  /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
     275  else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
     276  else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
     277   
     278  if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
     279 
     280  maxE=std::min(emax_model*1.001,maxE);
     281       
     282  G4double r = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
     283 
     284  if (IsIon) {
     285        G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
     286        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
     287  }     
     288 
     289  G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
     290 
     291  if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
     292 
     293 
     294
     295  x=r1-r;
     296  x=std::max(r1-r,0.001*mm);
     297 
     298  return x;
     299 
     300 
     301}
     302#include "G4EmCorrections.hh"
     303///////////////////////////////////////////////////////
     304//
     305
     306void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
     307{
     308
     309  G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
     310  if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
     311}

  • trunk/source/processes/electromagnetic/adjoint/src/G4InversePEEffect.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4InversePEEffect.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4InversePEEffect.hh"
    2730#include "G4VEmAdjointModel.hh"
    … …  
    3134//
    3235G4InversePEEffect::G4InversePEEffect(G4String process_name,G4AdjointPhotoElectricModel* aModel):
    33                                 G4VAdjointInverseScattering(process_name,false)
     36                                G4VAdjointReverseReaction(process_name,false)
    3437{theAdjointEMModel = aModel;
    35  theAdjointEMModel->SetSecondPartOfSameType(false);
     38 theAdjointEMModel->SetSecondPartOfSameType(false);
     39 SetIntegralMode(false);
     40 
     41 
    3642}
    3743/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4VEmAdjointModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4VEmAdjointModel.hh"
    2730#include "G4AdjointCSManager.hh"
    28 
    29 
    3031#include "G4Integrator.hh"
    3132#include "G4TrackStatus.hh"
    … …  
    3334#include "G4AdjointElectron.hh"
    3435#include "G4AdjointInterpolator.hh"
     36#include "G4PhysicsTable.hh"
    3537
    3638////////////////////////////////////////////////////////////////////////////////
    … …  
    3941name(nam)
    4042// lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
    41 { G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
    42   CorrectWeightMode =true;
    43   UseMatrix =true;
    44   UseMatrixPerElement = true;
    45   ApplyCutInRange = true;
    46   ApplyBiasing = true;
    47   UseOnlyOneMatrixForAllElements = true;
    48   IsIonisation =true;
    49   CS_biasing_factor =1.;
    50   //ApplyBiasing = false;
     43{
     44  G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
     45  second_part_of_same_type =false;
     46  theDirectEMModel=0;
    5147}
    5248////////////////////////////////////////////////////////////////////////////////
    … …  
    5450G4VEmAdjointModel::~G4VEmAdjointModel()
    5551{;}
    56 ////////////////////////////////////////////////////////////////////////////////
    57 //
    58 void G4VEmAdjointModel::SampleSecondaries(const G4Track& aTrack,
    59                        G4bool IsScatProjToProjCase,
    60                        G4ParticleChange* fParticleChange)
    61 {
    62 
    63   const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
    64   //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
    65   size_t ind=0;
    66   if (!UseMatrixPerElement) ind = currentMaterialIndex;
    67   //G4cout<<theAdjointPrimary<<std::endl;
    68   else if (!UseOnlyOneMatrixForAllElements) { //Select Material
    69         std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
    70         if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
    71         G4double rand_var= G4UniformRand();
    72         G4double SumCS=0.;
    73         for (size_t i=0;i<CS_Vs_Element->size();i++){
    74                 SumCS+=(*CS_Vs_Element)[i];
    75                 if (rand_var<=SumCS/lastCS){
    76                         ind=i;
    77                         break;
    78                 }
    79         }
    80         ind = currentMaterial->GetElement(ind)->GetIndex();
    81   }
    82  
    83  
    84  
    85  //Elastic inverse scattering //not correct in all the cases
    86  //---------------------------------------------------------
    87  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
    88  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
    89  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
    90  //G4cout<<adjointPrimKinEnergy<<std::endl;
    91  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
    92         return;
    93  }
    94  
    95  //Sample secondary energy
    96  //-----------------------
    97  G4double projectileKinEnergy;
    98 // if (!IsIonisation  ) {
    99         projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
    100                                                   adjointPrimKinEnergy,
    101                                                   IsScatProjToProjCase);
    102  //}
    103  /*else {
    104         projectileKinEnergy =   SampleAdjSecEnergyFromDiffCrossSectionPerAtom(adjointPrimKinEnergy,IsScatProjToProjCase);
    105         //G4cout<<projectileKinEnergy<<std::endl;
    106  }*/   
    107  //Weight correction
    108  //-----------------------
    109  
    110  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);                                     
    111  
    112                  
    113  
    114  
    115  
    116  //Kinematic
    117  //---------
    118  
    119  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
    120  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
    121  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
    122  
    123  
    124  
    125  //Companion
    126  //-----------
    127  G4double companionM0;
    128  companionM0=(adjointPrimTotalEnergy-adjointPrimKinEnergy);
    129  if (IsScatProjToProjCase) {
    130         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
    131  }
    132  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
    133  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;   
    134  
    135  
    136  //Projectile momentum
    137  //--------------------
    138  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
    139  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
    140  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
    141  G4double phi =G4UniformRand()*2.*3.1415926;
    142  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
    143  projectileMomentum.rotateUz(dir_parallel);
    144  
    145  
    146  
    147  if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
    148         fParticleChange->ProposeTrackStatus(fStopAndKill);
    149         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
    150         //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
    151  }
    152  else {
    153         fParticleChange->ProposeEnergy(projectileKinEnergy);
    154         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
    155  }     
    156 }
    157 ////////////////////////////////////////////////////////////////////////////////
    158 //
    159 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,  G4double , G4double )
    160 {
    161  G4double new_weight=old_weight;
    162  if (CorrectWeightMode) {
    163         G4double w_corr =1./CS_biasing_factor;
    164         //G4cout<<w_corr<<std::endl;
    165        
    166         /*G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(theAdjEquivOfDirectPrimPartDef,
    167                                                                                theAdjEquivOfDirectSecondPartDef,
    168                                                                                adjointPrimKinEnergy,projectileKinEnergy,
    169                                                                                aTrack.GetMaterialCutsCouple());
    170         w_corr = projectileKinEnergy;
    171         G4double Emin,Emax;
    172         if (IsScatProjToProjCase) {
    173                 Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
    174                 Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy, currentTcutForDirectSecond);
    175                
    176         }
    177         else {
    178                 Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
    179                 Emin = GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);   
    180         }
    181         w_corr *=std::log(Emax/Emin)/(Emax-Emin); */   
    182                
    183         new_weight*=w_corr;
    184  }
    185  G4cout<< "new weight"<<new_weight<<std::endl;
    186  fParticleChange->SetParentWeightByProcess(false);
    187  fParticleChange->SetSecondaryWeightByProcess(false);
    188  fParticleChange->ProposeParentWeight(new_weight);     
    189 }
    19052////////////////////////////////////////////////////////////////////////////////
    19153//                             
    … …  
    19557{
    19658  DefineCurrentMaterial(aCouple);
    197   //G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(G4AdjointElectron::AdjointElectron(),primEnergy,aCouple);
    198   //G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(G4AdjointElectron::AdjointElectron(), primEnergy,aCouple);
    199   if (IsScatProjToProjCase){
    200         lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
     59  preStepEnergy=primEnergy;
     60 
     61  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
     62  if (IsScatProjToProjCase)   CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
     63  lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
    20164                                                                        this,
    20265                                                                        primEnergy,
    20366                                                                        currentTcutForDirectSecond,
    204                                                                         true,
    205                                                                         CS_Vs_ElementForScatProjToProjCase);
    206         /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectPrimPartDef,primEnergy,aCouple);
    207         G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectPrimPartDef, primEnergy,aCouple);
    208         */
    209         //if (adjCS >0 )lastCS *=fwdCS/adjCS;
    210  
    211   }
    212   else {
    213         lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
    214                                                                         this,
    215                                                                         primEnergy,
    216                                                                         currentTcutForDirectSecond,
    217                                                                         false,
    218                                                                         CS_Vs_ElementForProdToProjCase);
    219         /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectSecondPartDef,primEnergy,aCouple);
    220         G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectSecondPartDef, primEnergy,aCouple);
    221         */
    222         //if (adjCS >0 )lastCS *=fwdCS/adjCS;                                                           
    223         //lastCS=0.;                                                           
    224   }
    225  
    226  
     67                                                                        IsScatProjToProjCase,
     68                                                                        *CS_Vs_Element);
     69  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
     70  else lastAdjointCSForProdToProjCase =lastCS;                                                         
     71       
     72 
     73 
    22774  return lastCS;
    22875                                                                       
    … …  
    23077////////////////////////////////////////////////////////////////////////////////
    23178//
    232 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
     79//General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
    23380G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
    23481                                      G4double kinEnergyProj,
    … …  
    24592        G4double Tmax=kinEnergyProj;
    24693        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
    247         return Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
    248         //it could be thta Tmax here should be DBLMAX
    249         //Tmax=DBLMAX;
    250        
    251         G4double E1=kinEnergyProd;
     94       
     95        G4double E1=kinEnergyProd;
    25296        G4double E2=kinEnergyProd*1.000001;
    25397        G4double dE=(E2-E1);
    … …  
    256100       
    257101        dSigmadEprod=(sigma1-sigma2)/dE;
    258         if (dSigmadEprod>1.) {
    259                 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<std::endl;
    260                 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<std::endl;
    261                 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<std::endl;
    262                
    263         }
    264        
    265        
    266        
    267102 }
    268  
    269  
    270 
    271        
    272103 return dSigmadEprod;   
    273104 
    … …  
    307138        G4double Tmax=kinEnergyProj;
    308139        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
    309         //it could be thta Tmax here should be DBLMAX
    310         //Tmax=DBLMAX;
    311        
    312         G4double E1=kinEnergyProd;
    313        
    314         G4double E2=kinEnergyProd*1.0001;
     140        G4double E1=kinEnergyProd;
     141        G4double E2=kinEnergyProd*1.0001;
    315142        G4double dE=(E2-E1);
    316143        G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
    317        
    318         //G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
    319         dSigmadEprod=sigma1/dE;
    320         if (dSigmadEprod <0) { //could happen with bremstrahlung dur to suppression effect
    321                 G4cout<<"Halllllllllllllllllllllllllllllllllllllllllllllllo "<<kinEnergyProj<<'\t'<<E1<<'\t'<<dSigmadEprod<<std::endl;
    322                 E1=kinEnergyProd;
    323                 E2=E1*1.1;
    324                 dE=E2-E1;
    325                 sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e50);
    326                 G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
    327                 dSigmadEprod=(sigma1-sigma2)/dE;
    328                 G4cout<<dSigmadEprod<<std::endl;
    329         }       
    330        
    331        
     144        G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
     145        dSigmadEprod=(sigma1-sigma2)/dE;
    332146 }
    333 
    334 
    335        
    336147 return dSigmadEprod;   
    337148 
    … …  
    356167//
    357168G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
    358   //return kinEnergyProj*kinEnergyProj;
    359   //ApplyBiasing=false;
     169 
     170 
    360171  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;   
    361   if (!ApplyBiasing) bias_factor =CS_biasing_factor;
    362   //G4cout<<bias_factor<<std::endl;
     172
     173
    363174  if (UseMatrixPerElement ) {
    364175        return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
    365176  }
    366   else {
     177  else  {
    367178        return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
    368179  }     
    369180}
    370 //////////////////////////////////////////////////////////////////////////////
    371 //
    372 G4double G4VEmAdjointModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
    373   G4double electron_mass_c2=0.51099906*MeV;
    374   G4double energy = kinEnergyProj + electron_mass_c2;
    375   G4double x   = kinEnergyProd/kinEnergyProj;
    376   G4double gam    = energy/electron_mass_c2;
    377   G4double gamma2 = gam*gam;
    378   G4double beta2  = 1.0 - 1.0/gamma2;
    379  
    380   G4double g = (2.0*gam - 1.0)/gamma2;
    381   G4double y = 1.0 - x;
    382   G4double fac=twopi_mc2_rcl2/electron_mass_c2;
    383   G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
    384   return dCS/kinEnergyProj;
    385  
    386  
    387 
    388 } 
     181
    389182////////////////////////////////////////////////////////////////////////////////
    390183//
    391184G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
    392   //return kinEnergyProj*kinEnergyProj;
    393   G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
    394   //ApplyBiasing=false;
    395  if (!ApplyBiasing) bias_factor = CS_biasing_factor;
    396  //G4cout<<bias_factor<<std::endl;
     185 
     186 G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
    397187 if (UseMatrixPerElement ) {
    398188        return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
    399189 }     
    400  else {
     190 else { 
    401191        return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
    402192 
    403193 }     
    404194               
     195}
     196////////////////////////////////////////////////////////////////////////////////
     197//
     198
     199G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
     200{
     201  return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
    405202}
    406203////////////////////////////////////////////////////////////////////////////////
    … …  
    411208                                G4double A ,
    412209                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
    413 { G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
    414   ASelectedNucleus= G4int(A);
    415   ZSelectedNucleus=G4int(Z);
     210{
     211  G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
     212  ASelectedNucleus= int(A);
     213  ZSelectedNucleus=int(Z);
    416214  kinEnergyProdForIntegration = kinEnergyProd;
    417215 
    … …  
    422220  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
    423221  G4double E1=minEProj;
    424   std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
    425   std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
     222  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
     223  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
    426224  log_ESec_vector->clear();
    427225  log_Prob_vector->clear();
    … …  
    429227  log_Prob_vector->push_back(-50.);
    430228 
    431   G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
     229  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
    432230  G4double fE=std::pow(10.,1./nbin_pro_decade);
    433231  G4double int_cross_section=0.;
    … …  
    436234 
    437235  while (E1 <maxEProj*0.9999999){
    438         //G4cout<<E1<<'\t'<<E2<<std::endl;
    439        
    440         int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
    441         //G4cout<<"int_cross_section 1 "<<'\t'<<int_cross_section<<std::endl;
     236        //G4cout<<E1<<'\t'<<E2<<G4endl;
     237       
     238        int_cross_section +=integral.Simpson(this,
     239        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
    442240        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
    443241        log_Prob_vector->push_back(std::log(int_cross_section));       
    … …  
    463261                                G4double A ,
    464262                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
    465 { G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
    466   ASelectedNucleus=G4int(A);
    467   ZSelectedNucleus=G4int(Z);
     263{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
     264  ASelectedNucleus=int(A);
     265  ZSelectedNucleus=int(Z);
    468266  kinEnergyScatProjForIntegration = kinEnergyScatProj;
    469267 
    … …  
    479277 
    480278 
    481   std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
    482   std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
     279  std::vector< double>*  log_ESec_vector = new std::vector< double>();
     280  std::vector< double>*  log_Prob_vector = new std::vector< double>();
    483281  log_ESec_vector->push_back(std::log(dEmin));
    484282  log_Prob_vector->push_back(-50.);
    485   G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
     283  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
    486284  G4double fE=std::pow(dEmax/dEmin,1./nbins);
     285 
     286 
     287 
     288 
    487289 
    488290  G4double int_cross_section=0.;
    … …  
    491293        dE2=dE1*fE;
    492294        int_cross_section +=integral.Simpson(this,
    493         &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
    494         //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
    495         log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
     295        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
     296        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
     297        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
    496298        log_Prob_vector->push_back(std::log(int_cross_section));       
    497299        dE1=dE2;
    498300
    499301  }
    500   /*G4cout<<"total int_cross_section"<<'\t'<<int_cross_section<<std::endl;
    501   G4cout<<"energy "<<kinEnergyScatProj<<std::endl;*/
    502  
    503  
    504302 
    505303 
    … …  
    519317                                G4double kinEnergyProd,
    520318                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
    521 { G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
     319{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
    522320  SelectedMaterial= aMaterial;
    523321  kinEnergyProdForIntegration = kinEnergyProd;
    524   //G4cout<<aMaterial->GetName()<<std::endl;
    525   //G4cout<<kinEnergyProd/MeV<<std::endl;
    526   //compute the vector of integrated cross sections
     322   //compute the vector of integrated cross sections
    527323  //-------------------
    528324 
    … …  
    530326  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
    531327  G4double E1=minEProj;
    532   std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
    533   std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
     328  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
     329  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
    534330  log_ESec_vector->clear();
    535331  log_Prob_vector->clear();
    … …  
    537333  log_Prob_vector->push_back(-50.);
    538334 
    539   G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
     335  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
    540336  G4double fE=std::pow(10.,1./nbin_pro_decade);
    541337  G4double int_cross_section=0.;
    … …  
    544340 
    545341  while (E1 <maxEProj*0.9999999){
    546         //G4cout<<E1<<'\t'<<E2<<std::endl;
    547        
    548         int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
    549         //G4cout<<"int_cross_section 1 "<<E1<<'\t'<<int_cross_section<<std::endl;
     342       
     343        int_cross_section +=integral.Simpson(this,
     344        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
    550345        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
    551346        log_Prob_vector->push_back(std::log(int_cross_section));
    … …  
    557352  res_mat.clear();
    558353 
    559   //if (int_cross_section >0.) {
     354  if (int_cross_section >0.) {
    560355        res_mat.push_back(log_ESec_vector);
    561356        res_mat.push_back(log_Prob_vector);
    562   //}   
     357  }
     358 
     359 
    563360 
    564361  return res_mat;
    … …  
    571368                                G4double kinEnergyScatProj,
    572369                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
    573 { G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
     370{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
    574371  SelectedMaterial= aMaterial;
    575372  kinEnergyScatProjForIntegration = kinEnergyScatProj;
    576   /*G4cout<<name<<std::endl;
    577   G4cout<<aMaterial->GetName()<<std::endl;
    578   G4cout<<kinEnergyScatProj/MeV<<std::endl;*/
     373 
    579374  //compute the vector of integrated cross sections
    580375  //-------------------
    … …  
    590385 
    591386 
    592   std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
    593   std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
     387  std::vector< double>*  log_ESec_vector = new std::vector< double>();
     388  std::vector< double>*  log_Prob_vector = new std::vector< double>();
    594389  log_ESec_vector->push_back(std::log(dEmin));
    595390  log_Prob_vector->push_back(-50.);
    596   G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
     391  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
    597392  G4double fE=std::pow(dEmax/dEmin,1./nbins);
    598393 
    … …  
    602397        dE2=dE1*fE;
    603398        int_cross_section +=integral.Simpson(this,
    604         &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
    605         //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
    606         log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
     399        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
     400        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
    607401        log_Prob_vector->push_back(std::log(int_cross_section));
    608402        dE1=dE2;
    … …  
    631425  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
    632426  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
    633   std::vector< G4double >* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
    634   //G4double dLog = theMatrix->GetDlog();
    635  
    636  
     427  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
    637428 
    638429  if (theLogPrimEnergyVector->size() ==0){
    639         G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
    640         G4cout<<"The sampling procedure will be stopped."<<std::endl;
     430        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
     431        G4cout<<"The sampling procedure will be stopped."<<G4endl;
    641432        return 0.;
    642433       
    … …  
    650441  G4double aLogPrimEnergy1,aLogPrimEnergy2;
    651442  G4double aLogCS1,aLogCS2;
    652    G4double log01,log02;
    653   std::vector< G4double>* aLogSecondEnergyVector1 =0;
    654   std::vector< G4double>* aLogSecondEnergyVector2  =0;
    655   std::vector< G4double>* aLogProbVector1=0;
    656   std::vector< G4double>* aLogProbVector2=0;
     443  G4double log01,log02;
     444  std::vector< double>* aLogSecondEnergyVector1 =0;
     445  std::vector< double>* aLogSecondEnergyVector2  =0;
     446  std::vector< double>* aLogProbVector1=0;
     447  std::vector< double>* aLogProbVector2=0;
    657448  std::vector< size_t>* aLogProbVectorIndex1=0;
    658449  std::vector< size_t>* aLogProbVectorIndex2=0;
    … …  
    674465  G4double Emax=0.;
    675466  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
    676         //G4cout<<"Here "<<std::endl;
    677         if (ApplyCutInRange) {
     467        Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
     468        Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
     469        G4double dE=0;
     470        if (Emin < Emax ){
     471            if (ApplyCutInRange) {
    678472                if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
    679                 /*if (IsIonisation){
    680                         G4double inv_Tcut= 1./currentTcutForDirectSecond;
    681                         G4double inv_dE=inv_Tcut-rand_var*(inv_Tcut-1./aPrimEnergy);
    682                         Esec= aPrimEnergy+1./inv_dE;
    683                         //return Esec;
    684                         G4double dE1=currentTcutForDirectSecond;
    685                         G4double dE2=currentTcutForDirectSecond*1.00001;
    686                         G4double dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
    687                         G4double dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
    688                         G4double alpha1=std::log(dCS1/dCS2)/std::log(dE1/dE2);
    689                         G4double a1=dCS1/std::pow(dE1,alpha1);
    690                         dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
    691                         dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
    692                        
    693                         return Esec;
    694                        
    695                        
    696                        
    697                         dE1=aPrimEnergy/1.00001;
    698                         dE2=aPrimEnergy;
    699                         dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
    700                         dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
    701                         G4double alpha2=std::log(dCS1/dCS2)/std::log(dE1/dE2);
    702                         G4double a2=dCS1/std::pow(dE1,alpha1);
    703                         return Esec;
    704                        
    705                        
    706                        
    707                        
    708                 }*/     
     473               
    709474                log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
    710475                log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
    711476               
    712         }       
    713         log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
    714         log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
    715        
    716         /*log_dE1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
    717         log_dE2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
    718         */                                 
    719                                            
    720        
    721        
    722        
    723        
    724         Esec = aPrimEnergy +
    725         std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
    726        
    727         Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy);
    728         Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
     477            }   
     478            log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
     479            log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
     480             dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
     481        }
     482       
     483        Esec = aPrimEnergy +dE;
    729484        Esec=std::max(Esec,Emin);
    730485        Esec=std::min(Esec,Emax);
    731486       
    732        
    733         //G4cout<<"Esec "<<Esec<<std::endl;
    734         //if (Esec > 2.*aPrimEnergy && second_part_of_same_type) Esec = 2.*aPrimEnergy;
    735487  }
    736488  else { //Tcut condition is already full-filled
    737         /*G4cout<<"Start "<<std::endl;
    738         G4cout<<std::exp((*aLogProbVector1)[0])<<std::endl;
    739         G4cout<<std::exp((*aLogProbVector2)[0])<<std::endl;*/
    740         /*G4double inv_E1= .5/aPrimEnergy;
    741                
    742         G4double inv_E=inv_E1-rand_var*(inv_E1-0.00001);
    743         Esec= 1./inv_E;
    744         return Esec;*/
     489       
    745490        log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
    746491        log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
    747         /*log_E1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
    748         log_E2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
    749         */
    750        
    751        
    752         /*G4cout<<std::exp(log_E1)<<std::endl;
    753         G4cout<<std::exp(log_E2)<<std::endl;*/
    754492       
    755493        Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
    … …  
    767505 
    768506                                                                                   
     507}
     508
     509//////////////////////////////////////////////////////////////////////////////
     510//                             
     511G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
     512{ SelectCSMatrix(IsScatProjToProjCase);
     513  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
     514}
     515//////////////////////////////////////////////////////////////////////////////
     516//
     517void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
     518{
     519  indexOfUsedCrossSectionMatrix=0;
     520  if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
     521  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
     522        std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
     523        lastCS=lastAdjointCSForScatProjToProjCase;
     524        if ( !IsScatProjToProjCase) {
     525                CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
     526                lastCS=lastAdjointCSForProdToProjCase;
     527        }       
     528        G4double rand_var= G4UniformRand();
     529        G4double SumCS=0.;
     530        size_t ind=0;
     531        for (size_t i=0;i<CS_Vs_Element->size();i++){
     532                SumCS+=(*CS_Vs_Element)[i];
     533                if (rand_var<=SumCS/lastCS){
     534                        ind=i;
     535                        break;
     536                }
     537        }
     538        indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
     539  }
    769540}
    770541//////////////////////////////////////////////////////////////////////////////
    … …  
    801572        do {
    802573                q = G4UniformRand();
    803                 x = std::pow(xmin, q);
     574                x = pow(xmin, q);
    804575                E=x*Emax;
    805576                greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
    … …  
    814585 
    815586  return E;
    816  
    817  
    818  
    819  
    820                                                                                    
     587}
     588
     589////////////////////////////////////////////////////////////////////////////////
     590//
     591void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
     592                                                  G4double old_weight, 
     593                                                  G4double adjointPrimKinEnergy,
     594                                                  G4double projectileKinEnergy,
     595                                                  G4bool IsScatProjToProjCase)
     596{
     597 G4double new_weight=old_weight;
     598 G4double w_corr =1./CS_biasing_factor;
     599 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
     600 
     601 
     602 lastCS=lastAdjointCSForScatProjToProjCase;
     603 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
     604 if (adjointPrimKinEnergy !=preStepEnergy){ //Is that in all cases needed???
     605        G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
     606                                                 ,IsScatProjToProjCase );
     607        w_corr*=post_stepCS/lastCS;                     
     608 }
     609       
     610 new_weight*=w_corr;
     611
     612 //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
     613 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
     614                                                        //by the factor adjointPrimKinEnergy/projectileKinEnergy
     615   
     616
     617
     618 fParticleChange->SetParentWeightByProcess(false);
     619 fParticleChange->SetSecondaryWeightByProcess(false);
     620 fParticleChange->ProposeParentWeight(new_weight);     
    821621}
    822622//////////////////////////////////////////////////////////////////////////////
    … …  
    830630//
    831631G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
    832 { return PrimAdjEnergy+Tcut;
     632{ G4double Emin=PrimAdjEnergy;
     633  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
     634  return Emin;
    833635}
    834636//////////////////////////////////////////////////////////////////////////////
    … …  
    853655        currentMaterialIndex = currentMaterial->GetIndex();
    854656        size_t idx=56;
    855    
     657        currentTcutForDirectPrim =0.00000000001;
    856658        if (theAdjEquivOfDirectPrimPartDef) {
    857659                if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
    858660                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
    859661                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
    860                 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
    861                 currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
     662                if (idx <56){
     663                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
     664                        currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
     665                }       
    862666        }       
    863667   
    864    
     668        currentTcutForDirectSecond =0.00000000001;
    865669        if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
    866670                currentTcutForDirectSecond = currentTcutForDirectPrim;
    … …  
    873677                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
    874678                        currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
     679                        if (idx <56){
     680                                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
     681                                currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
     682                        }
     683                       
     684                       
    875685                }
    876686        }
    877687  }
    878688}
     689////////////////////////////////////////////////////////////////////////////////////////////
     690//
     691void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
     692{ HighEnergyLimit=aVal;
     693  if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
     694}
     695////////////////////////////////////////////////////////////////////////////////////////////
     696//
     697void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
     698{
     699  LowEnergyLimit=aVal;
     700  if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
     701}
     702////////////////////////////////////////////////////////////////////////////////////////////
     703//
     704void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
     705{
     706  theAdjEquivOfDirectPrimPartDef=aPart;
     707  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
     708                                        theDirectPrimaryPartDef=G4Electron::Electron();
     709  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
     710                                        theDirectPrimaryPartDef=G4Gamma::Gamma();
     711}

  • trunk/source/processes/electromagnetic/adjoint/src/G4eInverseBremsstrahlung.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseBremsstrahlung.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629#include "G4eInverseBremsstrahlung.hh"
    2730#include "G4VEmAdjointModel.hh"
    … …  
    3134//
    3235G4eInverseBremsstrahlung::G4eInverseBremsstrahlung(G4bool whichScatCase,G4String process_name,G4AdjointBremsstrahlungModel* aBremAdjointModel):
    33                                 G4VAdjointInverseScattering(process_name,whichScatCase)
     36                                G4VAdjointReverseReaction(process_name,whichScatCase)
    3437{theAdjointEMModel = aBremAdjointModel;
    35  theAdjointEMModel->SetSecondPartOfSameType(false);
     38 theAdjointEMModel->SetSecondPartOfSameType(false);
     39 if (IsScatProjToProjCase) SetIntegralMode(true);
     40 else   SetIntegralMode(false);
    3641}
    3742/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • trunk/source/processes/electromagnetic/adjoint/src/G4eInverseCompton.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseCompton.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629///////////////////////////////////////////////////////
    2730// File name:     G4eInverseCompton
    … …  
    3942//
    4043G4eInverseCompton::G4eInverseCompton(G4bool whichScatCase,G4String process_name,G4AdjointComptonModel* aComptonAdjointModel):
    41                                 G4VAdjointInverseScattering(process_name,whichScatCase)
     44                                G4VAdjointReverseReaction(process_name,whichScatCase)
    4245{theAdjointEMModel = aComptonAdjointModel;
    43  theAdjointEMModel->SetSecondPartOfSameType(false);
     46 theAdjointEMModel->SetSecondPartOfSameType(false);
     47 SetIntegralMode(false);
     48 /*if (IsScatProjToProjCase) SetIntegralMode(false);
     49 else   SetIntegralMode(true); */
    4450}
    4551/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • trunk/source/processes/electromagnetic/adjoint/src/G4eInverseIonisation.cc

    r966 r1196  
    2424// ********************************************************************
    2525//
     26// $Id: G4eInverseIonisation.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
     28//
    2629///////////////////////////////////////////////////////
    2730// File name:     G4eInverseIonisation
    … …  
    3841//
    3942G4eInverseIonisation::G4eInverseIonisation(G4bool whichScatCase,G4String process_name,G4VEmAdjointModel* aEmAdjointModel):
    40                                 G4VAdjointInverseScattering(process_name,whichScatCase)
     43                                G4VAdjointReverseReaction(process_name,whichScatCase)
    4144{theAdjointEMModel = aEmAdjointModel;
    4245 theAdjointEMModel->SetSecondPartOfSameType(true);
     46 SetIntegralMode(true);
     47
    4348}
    4449/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.