Changeset 1196 for trunk/source/processes/electromagnetic/adjoint/src

Timestamp:

Nov 25, 2009, 5:13:58 PM (15 years ago)

Author:

garnier

Message:

update CVS release candidate geant4.9.3.01

Location:

trunk/source/processes/electromagnetic/adjoint/src

Files:

• G4AdjointAlongStepWeightCorrection.cc (modified) (6 diffs)
• G4AdjointBremsstrahlungModel.cc (modified) (4 diffs)
• G4AdjointCSManager.cc (modified) (41 diffs)
• G4AdjointCSMatrix.cc (modified) (8 diffs)
• G4AdjointComptonModel.cc (modified) (15 diffs)
• G4AdjointInterpolator.cc (modified) (9 diffs)
• G4AdjointPhotoElectricModel.cc (modified) (12 diffs)
• G4ContinuousGainOfEnergy.cc (modified) (10 diffs)
• G4InversePEEffect.cc (modified) (2 diffs)
• G4VEmAdjointModel.cc (modified) (34 diffs)
• G4eInverseBremsstrahlung.cc (modified) (2 diffs)
• G4eInverseCompton.cc (modified) (2 diffs)
• G4eInverseIonisation.cc (modified) (2 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointAlongStepWeightCorrection.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointAlongStepWeightCorrection.hh"
 #include "G4Step.hh"
@@ -30 +33 @@
 #include "G4AdjointCSManager.hh"
 
+#ifdef WIN32
+#include <G4float .h>
+#endif
 
 ///////////////////////////////////////////////////////
@@ -45 +51 @@
 {; 
 }
-
-
 ///////////////////////////////////////////////////////
 //
-
 void G4AdjointAlongStepWeightCorrection::PreparePhysicsTable(
      const G4ParticleDefinition& )
 {
 ; 
-
 }
-
 ///////////////////////////////////////////////////////
 //
@@ -63 +64 @@
 {;
 }
-
-
-
-
 ///////////////////////////////////////////////////////
 //
@@ -85 +82 @@
                                                                         preStepKinEnergy,Tkin, currentCouple,length);
         
+  
  
+  G4double new_weight=weight_correction*track.GetWeight();
   
-  G4double new_weight=weight_correction*track.GetWeight();
+  //if (weight_correction >2.) new_weight=1.e-300;
+  
+  
+  //The following test check for zero weight.
+  //This happens after weight correction of gamma for photo electric effect.
+  //When the new weight is 0 it will be later on consider as nan by G4.
+  //Therefore we do put a lower limit of 1.e-300. for new_weight 
+  //Correction by L.Desorgher on 15 July 2009 
+#ifdef WIN32
+  if (!!_isnan(new_weight) || new_weight==0){
+#else
+  if (std::isnan(new_weight) || new_weight==0){
+#endif
+                //G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
+                new_weight=1.e-300;
+  }
+  
+  //G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
   fParticleChange->SetParentWeightByProcess(false);
   fParticleChange->SetSecondaryWeightByProcess(false);
@@ -96 +112 @@
 
 }
+///////////////////////////////////////////////////////
+//
+G4double G4AdjointAlongStepWeightCorrection::GetContinuousStepLimit(const G4Track& track,
+                G4double , G4double , G4double& )
+{ 
+  G4double x = DBL_MAX;
+  DefineMaterial(track.GetMaterialCutsCouple());
+  preStepKinEnergy = track.GetKineticEnergy();
+  return x;
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointBremsstrahlungModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointBremsstrahlungModel.hh"
 #include "G4AdjointCSManager.hh"
@@ -30 +33 @@
 #include "G4ParticleChange.hh"
 #include "G4AdjointElectron.hh"
+#include "G4AdjointGamma.hh"
+#include "G4Electron.hh"
+
 #include "G4Timer.hh"
+//#include "G4PenelopeBremsstrahlungModel.hh"
+
 
 ////////////////////////////////////////////////////////////////////////////////
 //
 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
- G4VEmAdjointModel("AdjointBremModel"),
- probsup(1.0),
- MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
- LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
- theLPMflag(true)
-
-{ isElectron= true;
-  SetUseMatrix(true);
+ G4VEmAdjointModel("AdjointeBremModel"),
+  MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi)
+{ 
+  SetUseMatrix(false);
   SetUseMatrixPerElement(false);
+  
+  theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel");
+  theDirectEMModel=theDirectStdBremModel;
+ // theDirectPenelopeBremModel =0;
+        
   SetApplyCutInRange(true);
-  SetIsIonisation(false);
   highKinEnergy= 100.*TeV;
   lowKinEnergy = 1.0*keV;
   theTimer =new G4Timer();
   
-  theTimer->Start();
-  InitialiseParameters();
-  theTimer->Stop();
-  G4cout<<"Time elapsed in second for the initialidation of AdjointBrem "<<theTimer->GetRealElapsed()<<std::endl;
-  
-  ModeldCS="MODEL1";
+  theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
+  theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
+  theDirectPrimaryPartDef=G4Electron::Electron();
+  second_part_of_same_type=false;
+  
+  /*UsePenelopeModel=false;
+  if (UsePenelopeModel) {
+        G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
+        theEmModelManagerForFwdModels = new G4EmModelManager();
+        isPenelopeModelInitialised = false;
+        G4VEmFluctuationModel* f=0;
+        G4Region* r=0;
+        theDirectEMModel=thePenelopeModel;
+        theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
+  }
+  */    
+  
+
   
 }
+
+////////////////////////////////////////////////////////////////////////////////
+//
+void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
+                       G4bool IsScatProjToProjCase,
+                       G4ParticleChange* fParticleChange)
+{
+ if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
+
+ const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
+ DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
+ 
+ 
+ G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
+ G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
+ 
+ if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
+        return;
+ }
+  
+  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
+                                                        IsScatProjToProjCase);
+ //Weight correction
+ //-----------------------                                         
+ CorrectPostStepWeight(fParticleChange, 
+                       aTrack.GetWeight(), 
+                       adjointPrimKinEnergy,
+                       projectileKinEnergy,
+                       IsScatProjToProjCase);   
+ 
+ 
+ //Kinematic
+ //---------
+ G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
+ G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
+ G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
+ G4double projectileP = std::sqrt(projectileP2);
+ 
+ 
+ //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
+ //------------------------------------------------
+  G4double u;
+  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
+
+  if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
+     else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
+
+  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
+
+  G4double sint = std::sin(theta);
+  G4double cost = std::cos(theta);
+
+  G4double phi = twopi * G4UniformRand() ;
+  
+  G4ThreeVector projectileMomentum;
+  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
+  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
+        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
+        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
+        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
+        G4double sint1 =  std::sqrt(1.-cost1*cost1);
+        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
+  
+  }
+  
+  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
+ 
+ 
+ 
+  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
+        fParticleChange->ProposeTrackStatus(fStopAndKill);
+        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
+  }
+  else {
+        fParticleChange->ProposeEnergy(projectileKinEnergy);
+        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
+        
+  }     
+} 
+////////////////////////////////////////////////////////////////////////////////
+//
+void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
+                       G4bool IsScatProjToProjCase,
+                       G4ParticleChange* fParticleChange)
+{ 
+
+ const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
+ DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
+ 
+ 
+ G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
+ G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
+ 
+ if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
+        return;
+ }
+  
+ G4double projectileKinEnergy =0.;
+ G4double gammaEnergy=0.;
+ G4double diffCSUsed=0.; 
+ if (!IsScatProjToProjCase){
+        gammaEnergy=adjointPrimKinEnergy;
+        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
+        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
+        if (Emin>=Emax) return;
+        projectileKinEnergy=Emin*pow(Emax/Emin,G4UniformRand());
+        diffCSUsed=lastCZ/projectileKinEnergy;
+        
+ }
+ else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
+        G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
+        if (Emin>=Emax) return;
+        G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
+        G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
+        //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
+        projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*pow(f2,G4UniformRand()));
+        gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
+        diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
+        
+ }
+  
+  
+  
+                                                        
+ //Weight correction
+ //-----------------------
+ //First w_corr is set to the ratio between adjoint total CS and fwd total CS
+ G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
+
+ //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
+ //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.
+ //Basically any other differential CS   diffCS could be used here (example Penelope). 
+ 
+ G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
+ w_corr*=diffCS/diffCSUsed;
+           
+ G4double new_weight = aTrack.GetWeight()*w_corr;
+ fParticleChange->SetParentWeightByProcess(false);
+ fParticleChange->SetSecondaryWeightByProcess(false);
+ fParticleChange->ProposeParentWeight(new_weight);
+ 
+ //Kinematic
+ //---------
+ G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
+ G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
+ G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
+ G4double projectileP = std::sqrt(projectileP2);
+ 
+ 
+ //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
+ //------------------------------------------------
+  G4double u;
+  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
+
+  if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
+     else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
+
+  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
+
+  G4double sint = std::sin(theta);
+  G4double cost = std::cos(theta);
+
+  G4double phi = twopi * G4UniformRand() ;
+  
+  G4ThreeVector projectileMomentum;
+  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
+  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
+        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
+        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
+        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
+        G4double sint1 =  std::sqrt(1.-cost1*cost1);
+        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
+  
+  }
+  
+  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
+ 
+ 
+ 
+  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
+        fParticleChange->ProposeTrackStatus(fStopAndKill);
+        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
+  }
+  else {
+        fParticleChange->ProposeEnergy(projectileKinEnergy);
+        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
+        
+  }     
+} 
 ////////////////////////////////////////////////////////////////////////////////
 //
 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
 {;}
-////////////////////////////////////////////////////////////////////////////////
-//
-/*G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(
-                                      const G4Material* aMaterial,
-                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
-                                      G4double kinEnergyProd // kinetic energy of the secondary particle 
-                                      )
-{
- 
- static const G4double
-     ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
-     ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
-     ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
-
-  static const G4double
-     bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
-     bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
-     bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
-
-  static const G4double
-     al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
-     al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
-     al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
-
-  static const G4double
-     bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
-     bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
-     bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
-
-  static const G4double tlow = 1.*MeV;
-  
-  G4double dCrossEprod=0.;
-  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
-  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
- 
- 
- if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
-        
-        G4double cross = 0.0;
-        
-
-        
-        G4double E1=kinEnergyProd;
-        G4double E2=kinEnergyProd*1.000000001;
-        G4double dE=(E2-E1);
-
-        const G4ElementVector* theElementVector = aMaterial->GetElementVector();
-        const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
-        G4double dum=0.;
-  
-        for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
-
-                G4double fac=
-                
-                cross += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
-                      kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
-        
-                
-                
-        } 
-        dCrossEprod=(cross1-cross2)/dE; //first term
-        
-        //Now come the correction
-        //-----------------------
-        
-        //First compute fsig for E1
-        //-------------------------
-        
-        
-        G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
-        G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
-                                             *(aMaterial->GetElectronDensity());
-
-        G4double fsig = 0.;
-        G4int nmax = 100;
-        G4double vmin=std::log(E1);
-        G4double vmax=std::log(kinEnergyProj) ;
-        G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
-        G4double u,fac,c,v,dv,y ;
-        if(nn > 0) {
-
-                dv = (vmax-vmin)/nn ;
-                v  = vmin-dv ;
-                for(G4int n=0; n<=nn; n++) {
-
-                        v += dv;  
-                        u = std::exp(v);              
-                        fac = SupressionFunction(aMaterial, kinEnergyProj, u);
-                        y = u/kinEnergyProj;
-                        fac *= (4.-4.*y+3.*y*y)/3.;
-                        fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
-
-                        if ((n==0)||(n==nn)) c=0.5;
-                        else    c=1. ;
-
-                        fac  *= c;
-                        fsig += fac;
-                }
-                y = E1/kinEnergyProj ;
-                fsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
-
-        } 
-        else {
-                fsig = 1.;
-        }
-        if (fsig > 1.) fsig = 1.;
-        
-        dCrossEprod*=fsig;
-        //return dCrossEprod;
-        //Now we  compute dfsig 
-        //-------------------------
-        G4double dfsig = 0.;
-        nn=20;
-        vmax=std::log(E2) ;
-        dv = (vmax-vmin)/nn ;
-        v  = vmin-dv ;
-        for(G4int n=0; n<=nn; n++) {
-                v += dv;  
-                u = std::exp(v);              
-                fac = SupressionFunction(aMaterial, kinEnergyProj, u);
-                y = u/kinEnergyProj;
-                fac *= (4.-4.*y+3.*y*y)/3.;
-                fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
-
-                if ((n==0)||(n==nn)) c=0.5;
-                else    c=1. ;
-
-                fac  *= c;
-                dfsig += fac;
-        }
-        y = E1/kinEnergyProj;
-        dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
-        
-        dCrossEprod+=dfsig*cross1/dE;
-        
-        
-        
-         
-        
- }
- return dCrossEprod;
-  
-} 
-*/
+
+////////////////////////////////////////////////////////////////////////////////
+//
 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
                                       G4double kinEnergyProd // kinetic energy of the secondary particle 
                                       )
-{if (ModeldCS=="MODEL2") return  DiffCrossSectionPerVolumePrimToSecond2(aMaterial,
-                                                                 kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
+{/*if (UsePenelopeModel && !isPenelopeModelInitialised) {
+        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
+        isPenelopeModelInitialised =true;
+ } 
+ */                                                             
+ return  DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
+                                                                 kinEnergyProj, 
                                                                  kinEnergyProd);
- if (ModeldCS=="MODEL3") return  DiffCrossSectionPerVolumePrimToSecond3(aMaterial,
-                                                                 kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
-                                                                 kinEnergyProd);
- return  DiffCrossSectionPerVolumePrimToSecond1(aMaterial,
-                                                                 kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
-                                                                 kinEnergyProd);                                                                 
+ /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
+                                                                 kinEnergyProj, 
+                                                                 kinEnergyProd);*/                                                                                                                               
 }                                     
-////////////////////////////////////////////////////////////////////////////////
-// the one used till now
-G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond1(
+
+////////////////////////////////////////////////////////////////////////////////
+//
+G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
                                       const G4Material* aMaterial,
                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
@@ -233 +304 @@
  
  
+ //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
+ //This is what is applied in the discrete standard model before the  rejection test  that make a cooerction
+ //The application of the same rejection function is not possble here.
+ //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the 
+ // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and 
+ // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. 
+ // In the future we plan to use the brem secondary spectra from the G4Penelope implementation 
+ 
  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
-        
-        G4double cross1 = 0.0;
-        G4double cross2 = 0.0;
-
-        
-        G4double E1=kinEnergyProd;
-        G4double E2=kinEnergyProd*1.01;
-        G4double dE=(E2-E1);
-
-        const G4ElementVector* theElementVector = aMaterial->GetElementVector();
-        const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
-        G4double dum=0.;
-  
-        for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
-
-                cross1 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
-                      kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
-        
-                cross2 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
-                     kinEnergyProj, (*theElementVector)[i]->GetZ(), dum, E2);
-                
-        } 
-        dCrossEprod=(cross1-cross2)/dE; //first term
-        
-        //Now come the correction
-        //-----------------------
-        
-        //First compute fsig for E1
-        //-------------------------
-        
-        
-        G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
-        G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
-                                             *(aMaterial->GetElectronDensity());
-
-        G4double fsig1 = 0.;
-        G4int nmax = 100;
-        G4double vmin=std::log(E1);
-        G4double vmax=std::log(kinEnergyProj) ;
-        G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
-        G4double u,fac,c,v,dv,y ;
-        if(nn > 0) {
-
-                dv = (vmax-vmin)/nn ;
-                v  = vmin-dv ;
-                for(G4int n=0; n<=nn; n++) {
-
-                        v += dv;  
-                        u = std::exp(v);              
-                        fac = SupressionFunction(aMaterial, kinEnergyProj, u);
-                        y = u/kinEnergyProj;
-                        fac *= (4.-4.*y+3.*y*y)/3.;
-                        fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
-
-                        if ((n==0)||(n==nn)) c=0.5;
-                        else    c=1. ;
-
-                        fac  *= c;
-                        fsig1 += fac;
-                }
-                y = E1/kinEnergyProj ;
-                fsig1 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
-
-        } 
-        else {
-                fsig1 = 1.;
-        }
-        if (fsig1 > 1.) fsig1 = 1.;
-        
-        dCrossEprod*=fsig1;
-        
-        
-        G4double fsig2 = 0.;
-        vmin=std::log(E2);
-        nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
-        if(nn > 0) {
-
-                dv = (vmax-vmin)/nn ;
-                v  = vmin-dv ;
-                for(G4int n=0; n<=nn; n++) {
-
-                        v += dv;  
-                        u = std::exp(v);              
-                        fac = SupressionFunction(aMaterial, kinEnergyProj, u);
-                        y = u/kinEnergyProj;
-                        fac *= (4.-4.*y+3.*y*y)/3.;
-                        fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
-
-                        if ((n==0)||(n==nn)) c=0.5;
-                        else    c=1. ;
-
-                        fac  *= c;
-                        fsig2 += fac;
-                }
-                y = E2/kinEnergyProj ;
-                fsig2 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
-
-        } 
-        else {
-                fsig2 = 1.;
-        }
-        if (fsig2 > 1.) fsig2 = 1.;
-        
-
-        G4double dfsig=(fsig2-fsig1);
-        dCrossEprod+=dfsig*cross1/dE;
-        
-        dCrossEprod=(fsig1*cross1-fsig2*cross2)/dE;
-        
-        
-        
-        
-        
-        /*if (fsig < 1.){
-                //Now we  compute dfsig 
-                //-------------------------
-                G4double dfsig = 0.;
-                nn=20;
-                vmax=std::log(E2) ;
-                dv = (vmax-vmin)/nn ;
-                v  = vmin-dv ;
-                for(G4int n=0; n<=nn; n++) {
-                        v += dv;  
-                        u = std::exp(v);              
-                        fac = SupressionFunction(aMaterial, kinEnergyProj, u);
-                        y = u/kinEnergyProj;
-                        fac *= (4.-4.*y+3.*y*y)/3.;
-                        fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
-
-                        if ((n==0)||(n==nn)) c=0.5;
-                        else    c=1. ;
-
-                        fac  *= c;
-                        dfsig += fac;
-                }
-                y = E1/kinEnergyProj;
-                dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
-                dCrossEprod+=dfsig*cross1/dE;
-                
-        }       
-        */
-        
-        
-        
-        
-        
-         
-        
+        G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
+        dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
  }
  return dCrossEprod;
   
-} 
-
-
-////////////////////////////////////////////////////////////////////////////////
-//
-G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2(
-                                      const G4Material* aMaterial,
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
+                                      const G4Material* material,
                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
                                       G4double kinEnergyProd // kinetic energy of the secondary particle 
                                       )
 {
+ //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor 
+  //used in the direct model
+ 
  G4double dCrossEprod=0.;
- G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
- G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
- 
- 
- if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
-        
-        G4double dEdX1 = 0.0;
-        G4double dEdX2 = 0.0;
-
-        
-        G4double E1=kinEnergyProd;
-        G4double E2=kinEnergyProd*1.001;
-        G4double dE=(E2-E1);
-        //G4double dum=0.;
-        
-        dEdX1 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E1); 
-        dEdX2 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E2);
-        dCrossEprod=(dEdX2-dEdX1)/dE/E1;
-        
-         
-        
-        
-        
-         
-        
- }
+ 
+ const G4ElementVector* theElementVector = material->GetElementVector();
+ const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
+ G4double dum=0.;
+ G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
+ G4double dE=E2-E1;
+ for (size_t i=0; i<material->GetNumberOfElements(); i++) { 
+        G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
+        G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
+        dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
+   
+ }
+ 
+ //Now the Migdal correction
+ 
+ G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
+ G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
+                                             *(material->GetElectronDensity());
+                                             
+ 
+ G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
+                                                                    //model is different than the one used in the secondary sampling by a
+                                                                    //factor (1.+kp2) To be checked!
+ 
+ dCrossEprod*=MigdalFactor;
  return dCrossEprod;
   
@@ -428 +362 @@
 ////////////////////////////////////////////////////////////////////////////////
 //
-G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond3(
-                                      const G4Material* aMaterial,
-                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
-                                      G4double kinEnergyProd // kinetic energy of the secondary particle 
-                                      )
-{
- 
- return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
-                                                                 kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
-                                                                 kinEnergyProd);
-  
-}
-
-////////////////////////////////////////////////////////////////////////////////
-//
-G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material,
-                                 G4double kineticEnergy, G4double gammaEnergy)
-{
-  // supression due to the LPM effect+polarisation of the medium/
-  // supression due to the polarisation alone
-
-
-  G4double totEnergy = kineticEnergy+electron_mass_c2 ;
-  G4double totEnergySquare = totEnergy*totEnergy ;
-
-  G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
-
-  G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
-
-  G4double electronDensity = material->GetElectronDensity();
-
-  G4double sp = gammaEnergySquare/
-   (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
-
-  G4double supr = 1.0;
-
-  if (theLPMflag) {
-
-    G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
-
-    if (s2lpm < 1.) {
-
-      G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
-      G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
-      G4double splim = LPMgEnergyLimit2/
-        (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
-      G4double w = 1.+1./splim ;
-
-      if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
-      else                 w = s2lpm*(1.+1./sp);
-
-      supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ;
-      supr /= sp;    
-    } 
-    
-  } 
-  return supr;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-//
-void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
-                       G4bool IsScatProjToProjCase,
-                       G4ParticleChange* fParticleChange)
-{ 
-
-  //G4cout<<"Adjoint Brem"<<std::endl;
-  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
-  
-  size_t ind=0;
-  
-  if (UseMatrixPerElement ) { //Select Material
-        std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
-        if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
-        G4double rand_var= G4UniformRand();
-        G4double SumCS=0.;
-        for (size_t i=0;i<CS_Vs_Element->size();i++){
-                SumCS+=(*CS_Vs_Element)[i];
-                if (rand_var<=SumCS/lastCS){
-                        ind=i;
-                        break;
-                }
-        }
-  }
-  else  {
-        ind = currentMaterialIndex;
-  }
- 
- 
- //Elastic inverse scattering modified compared to general G4VEmAdjointModel
- //---------------------------
- G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
- G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
- //G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
- if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
-        return;
- }
- 
- //Sample secondary energy
- //-----------------------
- 
- G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
-                                                   adjointPrimKinEnergy,
-                                                   IsScatProjToProjCase);
-                                   
- 
- 
- 
- //Weight correction
- //-----------------------                                         
- CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);  
- 
- 
- //Kinematic
- //---------
- 
- G4double projectileM0 = electron_mass_c2;
- G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
- G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
- G4double projectileP = std::sqrt(projectileP2);
- 
- 
- //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
- //------------------------------------------------
-  G4double u;
-  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
-
-  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
-     else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
-
-  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
-
-  G4double sint = std::sin(theta);
-  G4double cost = std::cos(theta);
-
-  G4double phi = twopi * G4UniformRand() ;
-  
-  G4ThreeVector projectileMomentum;
-  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
-  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
-        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
-        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
-        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
-        G4double sint1 =  std::sqrt(1.-cost1*cost1);
-        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
-  
-  }
-  
-  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
- 
- 
- 
-  if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
-        fParticleChange->ProposeTrackStatus(fStopAndKill);
-        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
-        //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
+G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
+                                             G4double primEnergy,
+                                             G4bool IsScatProjToProjCase)
+{/* if (UsePenelopeModel && !isPenelopeModelInitialised) {
+        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
+        isPenelopeModelInitialised =true;
+  }
+  */
+  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
+  DefineCurrentMaterial(aCouple);
+  G4double Cross=0.;
+  lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/exp(1.));//this give the constant above
+  
+  if (!IsScatProjToProjCase ){
+        G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
+        G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
+        if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*log(Emax_proj/Emin_proj);
   }
   else {
-        fParticleChange->ProposeEnergy(projectileKinEnergy);
-        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
-        //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
-  }     
-} 
-////////////////////////////////////////////////////////////////////////////////
-//
-void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel)
-{theDirectBremModel=aModel;
- DefineDirectEMModel(aModel);
-} 
-////////////////////////////////////////////////////////////////////////////////
-//
-void G4AdjointBremsstrahlungModel::InitialiseParameters()
-{  
-  static const G4double
-     ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
-     ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
-     ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
-
-  static const G4double
-     bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
-     bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
-     bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
-
- /* static const G4double
-     al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
-     al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
-     al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
-
-  static const G4double
-     bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
-     bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
-     bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/
- 
-  
-  const G4ElementTable* theElementTable = G4Element::GetElementTable();
-  FZ.clear();
-  ah1.clear();
-  ah2.clear();
-  ah3.clear();
-  
-  bh1.clear();
-  bh2.clear();
-  bh3.clear();
-  
-  al0.clear();
-  al1.clear();
-  al2.clear();
-  
-  bl0.clear();
-  bl1.clear();
-  bl2.clear();
-  SigmaPerAtom.clear(); 
-  
-  for (size_t j=0; j<theElementTable->size();j++){
-        
-        G4Element* anElement=(*theElementTable)[j]; 
-        G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
-        FZ.push_back(lnZ* (4.- 0.55*lnZ));
-        G4double ZZ = anElement->GetIonisation()->GetZZ3();
-        
-        ah1.push_back(ah10 + ZZ* (ah11 + ZZ* ah12));
-        ah2.push_back(ah20 + ZZ* (ah21 + ZZ* ah22));
-        ah3.push_back(ah30 + ZZ* (ah31 + ZZ* ah32));
-
-        bh1.push_back(bh10 + ZZ* (bh11 + ZZ* bh12));
-        bh2.push_back(bh20 + ZZ* (bh21 + ZZ* bh22));
-        bh3.push_back(bh30 + ZZ* (bh31 + ZZ* bh32));
-        /*SigmaPerAtom.push_back(theDirectEMModel->ComputeCrossSectionPerAtom(
-                                        theDirectPrimaryPartDef,GetHighEnergyLimit()/2., 
-                                        anElement->GetZ(),1.,GetLowEnergyLimit(),1.e20));*/
-        
-        
+        G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
+        G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
+        if (Emax_proj>Emin_proj) Cross= lastCZ*log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
         
-  }     
-}
+  }
+  return Cross; 
+}                                            
+
+
+
+
+

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointCSManager.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointCSManager.hh"
 #include "G4AdjointCSMatrix.hh"
@@ -40 +43 @@
 #include "G4Electron.hh"
 #include "G4Gamma.hh"
+#include "G4Proton.hh"
 #include "G4AdjointElectron.hh"
 #include "G4AdjointGamma.hh"
+#include "G4AdjointProton.hh"
 #include "G4ProductionCutsTable.hh"
 #include "G4ProductionCutsTable.hh"
+#include <fstream>
+#include <iomanip>
 
 
@@ -65 +72 @@
   listOfForwardEmProcess.clear();
   listOfForwardEnergyLossProcess.clear();
-  theListOfAdjointParticlesInAction.clear();
+  theListOfAdjointParticlesInAction.clear(); 
+  EminForFwdSigmaTables.clear();
+  EminForAdjSigmaTables.clear();
+  EkinofFwdSigmaMax.clear();
+  EkinofAdjSigmaMax.clear();
   Tmin=0.1*keV;
   Tmax=100.*TeV;
-  nbins=240;
+  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)
   
   RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
   RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
+  RegisterAdjointParticle(G4AdjointProton::AdjointProton());
   
   verbose  = 1;
-  
-  consider_continuous_weight_correction =true;
-  consider_poststep_weight_correction =false;
+ 
+  lastPartDefForCS =0;
+  LastEkinForCS =0;
+  LastCSCorrectionFactor =1.;
+  
+  forward_CS_mode = true;
+  
+  currentParticleDef = 0;
+  
+  theAdjIon = 0;
+  theFwdIon = 0;  
+ 
  
 }
@@ -96 +117 @@
   if (anAdjPartDef && aProcess){
         RegisterAdjointParticle(anAdjPartDef);
-        int index=-1;
+        G4int index=-1;
         
         for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
@@ -111 +132 @@
   if (anAdjPartDef && aProcess){
         RegisterAdjointParticle(anAdjPartDef);
-        int index=-1;
+        G4int index=-1;
         for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
                 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
@@ -121 +142 @@
 //
 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
-{  int index=-1;
+{  G4int index=-1;
    for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
         if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
@@ -132 +153 @@
         listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
         theListOfAdjointParticlesInAction.push_back(aPartDef);
+        EminForFwdSigmaTables.push_back(std::vector<G4double> ());
+        EminForAdjSigmaTables.push_back(std::vector<G4double> ());
+        EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
+        EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
+        
    }
 }
@@ -162 +188 @@
         const G4ElementTable* theElementTable = G4Element::GetElementTable();
         const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
+        
+        G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
         for (size_t i=0; i<listOfAdjointEMModel.size();i++){
                 G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
-                G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<std::endl;
+                G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
                 if (aModel->GetUseMatrix()){
                         std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
@@ -173 +201 @@
                                 if (aModel->GetUseOnlyOneMatrixForAllElements()){
                                                 std::vector<G4AdjointCSMatrix*>
-                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 10);
+                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
                                                 aListOfMat1->push_back(two_matrices[0]);
                                                 aListOfMat2->push_back(two_matrices[1]);
@@ -180 +208 @@
                                         for (size_t j=0; j<theElementTable->size();j++){
                                                 G4Element* anElement=(*theElementTable)[j];
-                                                G4int Z = G4int(anElement->GetZ());
-                                                G4int A = G4int(anElement->GetA());
+                                                G4int Z = int(anElement->GetZ());
+                                                G4int A = int(anElement->GetA());
                                                 std::vector<G4AdjointCSMatrix*>
-                                                        two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 10);
+                                                        two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
                                                 aListOfMat1->push_back(two_matrices[0]);
                                                 aListOfMat2->push_back(two_matrices[1]);
@@ -193 +221 @@
                                         G4Material* aMaterial=(*theMaterialTable)[j];
                                         std::vector<G4AdjointCSMatrix*>
-                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 10);
+                                                two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
                                         aListOfMat1->push_back(two_matrices[0]);
                                         aListOfMat2->push_back(two_matrices[1]);
@@ -203 +231 @@
                         aModel->SetCSMatrices(aListOfMat1, aListOfMat2);        
                 }
-                else {  std::vector<G4AdjointCSMatrix*> two_empty_matrices;
+                else {  G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
+                        std::vector<G4AdjointCSMatrix*> two_empty_matrices;
                         theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
                         theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
@@ -209 +238 @@
                 }               
         }
-        G4cout<<"All adjoint cross section matrices are built "<<std::endl;
+        G4cout<<"              All adjoint cross section matrices are computed!"<<G4endl;
+        G4cout<<"======================================================================"<<G4endl;
+        
         CrossSectionMatrixesAreBuilt = true;
+
+
 }
 
@@ -221 +254 @@
   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
         G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
+        DefineCurrentParticle(thePartDef);
         theTotalForwardSigmaTableVector[i]->clearAndDestroy();
         theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
+        EminForFwdSigmaTables[i].clear();
+        EminForAdjSigmaTables[i].clear();
+        EkinofFwdSigmaMax[i].clear();
+        EkinofAdjSigmaMax[i].clear();
+        //G4cout<<thePartDef->GetParticleName();
+        
         for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
                 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
                 
+                /*
+                G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
+                G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
+                
+                std::fstream FileOutputAdjCS(file_name1, std::ios::out);
+                std::fstream FileOutputFwdCS(file_name2, std::ios::out);
+                
+                
+                
+                FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
+                FileOutputAdjCS<<std::setprecision(6);
+                FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
+                FileOutputFwdCS<<std::setprecision(6);  
+                */
+                          
+
                 //make first the total fwd CS table for FwdProcess
                 G4PhysicsVector* aVector =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
+                G4bool Emin_found=false;
+                size_t ind=0;
+                G4double sigma_max =0.;
+                G4double e_sigma_max =0.;
                 for(size_t l=0; l<aVector->GetVectorLength(); l++) { 
-                        G4double totCS=0;
+                        G4double totCS=0.;
                         G4double e=aVector->GetLowEdgeEnergy(l);
                         for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
@@ -235 +295 @@
                         }
                         for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
-                                totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e, couple);
-                        }
-                        //G4cout<<totCS<<std::endl;
+                                if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
+                                        size_t mat_index = couple->GetIndex();
+                                        G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
+                                        G4double chargeSqRatio =  currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
+                                        (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
+                                }
+                                G4double e1=e/massRatio;
+                                totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
+                        }
                         aVector->PutValue(l,totCS);
-                
-                }
+                        if (totCS>sigma_max){
+                                sigma_max=totCS;
+                                e_sigma_max = e;
+                                
+                        }
+                        //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
+                        
+                        if (totCS>0 && !Emin_found) {
+                                EminForFwdSigmaTables[i].push_back(e);
+                                Emin_found=true;
+                        }
+                        
+                
+                }
+                //FileOutputFwdCS.close();
+                
+                EkinofFwdSigmaMax[i].push_back(e_sigma_max);
+                
+                
+                if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
+                
                 theTotalForwardSigmaTableVector[i]->push_back(aVector);
                 
+                
+                Emin_found=false;
+                sigma_max=0;
+                e_sigma_max =0.;
+                ind=0;
                 G4PhysicsVector* aVector1 =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
                 for(size_t l=0; l<aVector->GetVectorLength(); l++) { 
                         G4double e=aVector->GetLowEdgeEnergy(l);
-                        G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e);
-                        //G4cout<<totCS<<std::endl;
+                        G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
                         aVector1->PutValue(l,totCS);
-                
-                }       
+                        if (totCS>sigma_max){
+                                sigma_max=totCS;
+                                e_sigma_max = e;
+                                
+                        }
+                        //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
+                        if (totCS>0 && !Emin_found) {
+                                EminForAdjSigmaTables[i].push_back(e);
+                                Emin_found=true;
+                        }
+                
+                }
+                //FileOutputAdjCS.close();
+                EkinofAdjSigmaMax[i].push_back(e_sigma_max);
+                if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
+                        
                 theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
                 
@@ -262 +365 @@
                                      const G4MaterialCutsCouple* aCouple)
 { DefineCurrentMaterial(aCouple);
-  int index=-1;
-  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
-        if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
-  }     
-  if (index == -1) return 0.;
-  
+  DefineCurrentParticle(aPartDef);      
   G4bool b;
-  return (((*theTotalAdjointSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
+  return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
   
   
@@ -279 +377 @@
                                      const G4MaterialCutsCouple* aCouple)
 { DefineCurrentMaterial(aCouple);
-  int index=-1;
-  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
-        if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
-  }     
-  if (index == -1) return 0.;
+  DefineCurrentParticle(aPartDef);
   G4bool b;
-  return (((*theTotalForwardSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
-  
-  
+  return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
+  
+  
+}
+                                     
+///////////////////////////////////////////////////////
+//                                   
+void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
+                                     const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
+{ DefineCurrentMaterial(aCouple);
+  DefineCurrentParticle(aPartDef);
+  emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
+  emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
+  
+  
+  
+}
+///////////////////////////////////////////////////////
+//                                   
+void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
+                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
+{ DefineCurrentMaterial(aCouple);
+  DefineCurrentParticle(aPartDef);
+  e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
+  G4bool b;
+  sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
+  e_sigma_max/=massRatio;
+  
+  
+}
+///////////////////////////////////////////////////////
+//                                   
+void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
+                                     const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
+{ DefineCurrentMaterial(aCouple);
+  DefineCurrentParticle(aPartDef);
+  e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
+  G4bool b;
+  sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
+  e_sigma_max/=massRatio;
+  
+  
+}                                    
+///////////////////////////////////////////////////////
+//                                   
+G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
+                                                                                G4double& fwd_TotCS)
+{ G4double corr_fac = 1.;
+  if (forward_CS_mode) {
+        fwd_TotCS=PrefwdCS;
+        if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
+                DefineCurrentMaterial(aCouple);
+                PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
+                PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
+                LastEkinForCS = PreStepEkin;
+                lastPartDefForCS = aPartDef;
+                if (PrefwdCS >0. &&  PreadjCS >0.) {
+                        forward_CS_is_used = true;
+                        LastCSCorrectionFactor = PrefwdCS/PreadjCS;
+                }
+                else {
+                        forward_CS_is_used = false;
+                        LastCSCorrectionFactor = 1.;
+                        
+                }
+                
+        }
+        corr_fac =LastCSCorrectionFactor;
+        
+        
+        
+  }  
+  else {
+        forward_CS_is_used = false;
+        LastCSCorrectionFactor = 1.;
+  }
+  fwd_TotCS=PrefwdCS;   
+  fwd_is_used = forward_CS_is_used;
+  return  corr_fac;
 }                                    
 ///////////////////////////////////////////////////////
@@ -293 +463 @@
 G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
                                      const G4MaterialCutsCouple* aCouple, G4double step_length)
-{ //G4double fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
-  
-  G4double corr_fac = 1.;
-  if (consider_continuous_weight_correction) {
-        
-        G4double adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
-        G4double PrefwdCS;
-        PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
-        G4double fwdCS = GetTotalForwardCS(aPartDef, (AfterStepEkin+PreStepEkin)/2.,aCouple);
-        G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
-        //if (aPartDef ==G4AdjointGamma::AdjointGamma()) G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
-        /*if (adjCS >0 ) corr_fac = std::exp((PrefwdCS-fwdCS)*step_length);
-        else corr_fac = std::exp(-fwdCS*step_length);*/
-        corr_fac *=std::exp((adjCS-fwdCS)*step_length);
-        corr_fac=std::max(corr_fac,1.e-6);
-        corr_fac *=PreStepEkin/AfterStepEkin;
-        
-  }
-  G4cout<<"Cont "<<corr_fac<<std::endl;
-  G4cout<<"Ekin0 "<<PreStepEkin<<std::endl;
-  G4cout<<"Ekin1 "<<AfterStepEkin<<std::endl;
-  G4cout<<"step_length "<<step_length<<std::endl;
+{  G4double corr_fac = 1.;
+  //return corr_fac;
+  //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
+  G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
+  G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
+  if (!forward_CS_is_used || pre_adjCS ==0. ||  after_fwdCS==0.) {
+        forward_CS_is_used=false;
+        G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
+        corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
+        LastCSCorrectionFactor = 1.;
+  }
+  else {
+        LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
+  }     
+        
+
+ 
   return corr_fac; 
 }                                    
 ///////////////////////////////////////////////////////
 //                                                           
-G4double G4AdjointCSManager::GetPostStepWeightCorrection(G4ParticleDefinition* , G4ParticleDefinition* ,
-                                                          G4double ,G4double ,
-                                                          const G4MaterialCutsCouple* )
-{ G4double corr_fac = 1.;
-  if (consider_poststep_weight_correction) {
-        /*G4double fwdCS = GetTotalForwardCS(aSecondPartDef, EkinPrim,aCouple);
-        G4double adjCS = GetTotalAdjointCS(aPrimPartDef, EkinPrim,aCouple);*/
-        //G4double fwd1CS = GetTotalForwardCS(aPrimPartDef, EkinPrim,aCouple);
-        //if (adjCS>0 && fwd1CS>0) adjCS = fwd1CS;
-        //corr_fac =fwdCS*EkinSecond/adjCS/EkinPrim;
-        //corr_fac = adjCS/fwdCS;
-  }
-  return corr_fac;
+G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
+{//return 1.; 
+ return  1./LastCSCorrectionFactor;
+        
 }                                                         
 ///////////////////////////////////////////////////////
 //
-double  G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
+G4double  G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
                                                          G4VEmAdjointModel* aModel, 
                                                          G4double PrimEnergy,
                                                          G4double Tcut,
                                                          G4bool IsScatProjToProjCase,
-                                                         std::vector<double>& CS_Vs_Element)
+                                                         std::vector<G4double>& CS_Vs_Element)
 { 
   
+  G4double EminSec=0;
+  G4double EmaxSec=0;
+  
+  if (IsScatProjToProjCase){
+        EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
+        EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
+  }
+  else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
+        EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
+        EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
+  }
+  if (EminSec >= EmaxSec) return 0.;
+
+
   G4bool need_to_compute=false;
   if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
@@ -381 +551 @@
         CS_Vs_Element.clear();
         if (!aModel->GetUseMatrix()){
-                return aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase);
+                CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
                                          
         
@@ -397 +567 @@
                                         CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
                         G4double factor=0.;
-                        for (size_t i=0;i<n_el;i++){
-                                size_t ind_el = aMaterial->GetElement(i)->GetIndex();
+                        for (size_t i=0;i<n_el;i++){ //this could be computed only once
+                                //size_t ind_el = aMaterial->GetElement(i)->GetIndex();
                                 factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
-                                G4AdjointCSMatrix* theCSMatrix;
-                                if (IsScatProjToProjCase){
-                                        theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
-                                }
-                                else    theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
-                                //G4double CS =0.;
-                                
-                                //G4cout<<CS<<std::endl;                        
-                                
                         }
                         CS *=factor;
@@ -417 +578 @@
                         for (size_t i=0;i<n_el;i++){
                                 size_t ind_el = aMaterial->GetElement(i)->GetIndex();
-                                //G4cout<<aMaterial->GetName()<<std::endl;
+                                //G4cout<<aMaterial->GetName()<<G4endl;
                                 G4AdjointCSMatrix* theCSMatrix;
                                 if (IsScatProjToProjCase){
@@ -426 +587 @@
                                 if (PrimEnergy > Tlow)
                                         CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
-                                //G4cout<<CS<<std::endl;                        
+                                //G4cout<<CS<<G4endl;                   
                                 CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i])); 
                         }
@@ -453 +614 @@
   G4double CS=0;
   for (size_t i=0;i<CS_Vs_Element.size();i++){
-        CS+=CS_Vs_Element[i];
-  }
-
+        CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
+        
+  }
   return CS;
-        
-  
-  
-  
-  
-  
-  
-  
 }                                                       
 ///////////////////////////////////////////////////////
@@ -473 +626 @@
                                                                    G4double Tcut,
                                                                    G4bool IsScatProjToProjCase)
-{ std::vector<double> CS_Vs_Element;
+{ std::vector<G4double> CS_Vs_Element;
   G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
   G4double rand_var= G4UniformRand();
@@ -498 +651 @@
 {
  G4double TotalCS=0.;
-// G4ParticleDefinition* theDirPartDef = GetForwardParticleEquivalent(aPartDef);
+
  DefineCurrentMaterial(aCouple);
-/* size_t idx=-1;
- if (theDirPartDef->GetParticleName() == "gamma") idx = 0;
- else if (theDirPartDef->GetParticleName() == "e-") idx = 1;
- else if (theDirPartDef->GetParticleName() == "e+") idx = 2;
- 
- //THe tCut computation is wrong this should be on Tcut per model the secondary determioming the Tcut
- const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
- //G4cout<<aVec<<std::endl;
- G4double Tcut =(*aVec)[aCouple->GetIndex()];*/
- //G4cout<<"Tcut "<<Tcut<<std::endl;
- //G4cout<<(*aVec)[0]<<std::endl;
-// G4double Tcut =converters[idx]->Convert(Rcut,aCouple->GetMaterial());
- 
-  
- std::vector<double> CS_Vs_Element;
+
+  
+ std::vector<G4double> CS_Vs_Element;
  for (size_t i=0; i<listOfAdjointEMModel.size();i++){
-        /*G4ParticleDefinition* theDirSecondPartDef = 
-                        GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
         
-        */
-        
-        
         G4double Tlow=0;
         if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
@@ -527 +663 @@
                 G4ParticleDefinition* theDirSecondPartDef = 
                         GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
-                G4int idx=-1;
+                size_t idx=56;
                 if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
                 else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
                 else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
-                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
-                Tlow =(*aVec)[aCouple->GetIndex()];
+                if (idx <56) {
+                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
+                        Tlow =(*aVec)[aCouple->GetIndex()];
+                }       
                 
         
@@ -539 +677 @@
         if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
                 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
-                        //G4cout<<"Yes1 before "<<std::endl;
                         TotalCS += ComputeAdjointCS(currentMaterial,
                                                        listOfAdjointEMModel[i], 
-                                                       Ekin, Tlow,true,CS_Vs_Element);
-                        //G4cout<<"Yes1 "<<Ekin<<'\t'<<TotalCS<<std::endl;                             
+                                                       Ekin, Tlow,true,CS_Vs_Element);                         
                 }
                 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
@@ -549 +685 @@
                                                        listOfAdjointEMModel[i], 
                                                        Ekin, Tlow,false, CS_Vs_Element);
-                                                       
-                        //G4cout<<"Yes2 "<<TotalCS<<std::endl;
                 }
                 
@@ -563 +697 @@
 std::vector<G4AdjointCSMatrix*>
 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
-                                                                        int nbin_pro_decade)
+                                                                        G4int nbin_pro_decade)
 { 
   G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
@@ -577 +711 @@
    
     
-   
-   
-   
-   
-   
    //Product to projectile backward scattering
    //-----------------------------------------
    G4double fE=std::pow(10.,1./nbin_pro_decade);
-   G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
+   G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    G4double E1=EkinMin;
    while (E1 <EkinMaxForProd){
         E1=std::max(EkinMin,E2);
         E1=std::min(EkinMaxForProd,E1);
-        std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
+        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
         if (aMat.size()>=2) {
-                std::vector< G4double >* log_ESecVec=aMat[0];
-                std::vector< G4double >* log_CSVec=aMat[1];
+                std::vector< double>* log_ESecVec=aMat[0];
+                std::vector< double>* log_CSVec=aMat[1];
                 G4double log_adjointCS=log_CSVec->back();
                 //normalise CSVec such that it becomes a probability vector
-                /*for (size_t j=0;j<log_CSVec->size();j++) (*log_CSVec)[j]=(*log_CSVec)[j]-log_adjointCS;
-                (*log_CSVec)[0]=-90.;*/
-        
-        
-                for (size_t j=0;j<log_CSVec->size();j++) {
-                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
-                        if (j==0) (*log_CSVec)[j] = 0.; 
-                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
-                        //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
+                for (size_t j=0;j<log_CSVec->size();j++) {
+                        if (j==0) (*log_CSVec)[j] = 0.; 
+                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
                 }       
-                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
+                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
                 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
         }       
@@ -616 +739 @@
    //-----------------------------------------
    
-   E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
+   E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    E1=EkinMin;
    while (E1 <EkinMaxForScat){
         E1=std::max(EkinMin,E2);
         E1=std::min(EkinMaxForScat,E1);
-        std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
+        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
         if (aMat.size()>=2) {
-                std::vector< G4double >* log_ESecVec=aMat[0];
-                std::vector< G4double >* log_CSVec=aMat[1];
+                std::vector< double>* log_ESecVec=aMat[0];
+                std::vector< double>* log_CSVec=aMat[1];
                 G4double log_adjointCS=log_CSVec->back();
                 //normalise CSVec such that it becomes a probability vector
                 for (size_t j=0;j<log_CSVec->size();j++) {
-                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
-                        if (j==0) (*log_CSVec)[j] = 0.; 
-                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
-                //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
+                        if (j==0) (*log_CSVec)[j] = 0.; 
+                        else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
                 }       
-                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
+                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
                 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
         }
@@ -641 +762 @@
    
    
-   
-   
-   
-   
-   
   std::vector<G4AdjointCSMatrix*> res;
   res.clear();
-  
   res.push_back(theCSMatForProdToProjBackwardScattering);
   res.push_back(theCSMatForScatProjToProjBackwardScattering);
   
 
-#ifdef TEST_MODE
+/*
   G4String file_name;
   std::stringstream astream;
@@ -662 +777 @@
   theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
  
-  /*G4AdjointCSMatrix* aMat1 = new G4AdjointCSMatrix(false);
-  G4AdjointCSMatrix* aMat2 = new G4AdjointCSMatrix(true);
-  
-  aMat1->Read(G4String("test_Z")+str_Z+"_1.txt");
-  aMat2->Read(G4String("test_Z")+str_Z+"_2.txt");
-  aMat1->Write(G4String("test_Z")+str_Z+"_11.txt");
-  aMat2->Write(G4String("test_Z")+str_Z+"_22.txt"); */
-#endif 
+*/
+
   
   return res;
@@ -702 +811 @@
    //-----------------------------------------
    G4double fE=std::pow(10.,1./nbin_pro_decade);
-   G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
+   G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    G4double E1=EkinMin;
    while (E1 <EkinMaxForProd){
         E1=std::max(EkinMin,E2);
         E1=std::min(EkinMaxForProd,E1);
-        std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
+        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
         if (aMat.size()>=2) {
-                std::vector< G4double >* log_ESecVec=aMat[0];
-                std::vector< G4double >* log_CSVec=aMat[1];
+                std::vector< double>* log_ESecVec=aMat[0];
+                std::vector< double>* log_CSVec=aMat[1];
                 G4double log_adjointCS=log_CSVec->back();
         
                 //normalise CSVec such that it becomes a probability vector
                 for (size_t j=0;j<log_CSVec->size();j++) {
-                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
+                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
                         if (j==0) (*log_CSVec)[j] = 0.; 
                         else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
-                        //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
+                        //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
                 }       
-                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
+                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
                 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
         }       
@@ -733 +842 @@
    //-----------------------------------------
    
-   E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
+   E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
    E1=EkinMin;
    while (E1 <EkinMaxForScat){
         E1=std::max(EkinMin,E2);
         E1=std::min(EkinMaxForScat,E1);
-        std::vector< std::vector< G4double >* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
+        std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
         if (aMat.size()>=2) {
-                std::vector< G4double >* log_ESecVec=aMat[0];
-                std::vector< G4double >* log_CSVec=aMat[1];
+                std::vector< double>* log_ESecVec=aMat[0];
+                std::vector< double>* log_CSVec=aMat[1];
                 G4double log_adjointCS=log_CSVec->back();
         
                 for (size_t j=0;j<log_CSVec->size();j++) {
-                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
+                        //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
                         if (j==0) (*log_CSVec)[j] = 0.; 
                         else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
-                //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
+                //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
+ 
                 }       
-                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
+                (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
         
                 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
@@ -770 +880 @@
   res.push_back(theCSMatForScatProjToProjBackwardScattering); 
   
-#ifdef TEST_MODE
+ /*
   theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
   theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
-#endif
+*/
 
 
@@ -786 +896 @@
 {
  if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
- if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
+ else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
+ else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
+ else if (theFwdPartDef ==theFwdIon) return theAdjIon;
+ 
  return 0;      
 }
@@ -794 +907 @@
 {
  if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
- if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
+ else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
+ else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
+ else if (theAdjPartDef == theAdjIon) return theFwdIon;
  return 0;      
 }
@@ -805 +920 @@
     currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
     currentMatIndex = couple->GetIndex();
-    //G4cout<<"Index material "<<currentMatIndex<<std::endl;
+    lastPartDefForCS =0;
+    LastEkinForCS =0;
+    LastCSCorrectionFactor =1.;
   }  
 }
 
-
-
-///////////////////////////////////////////////////////
-//
-double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
+///////////////////////////////////////////////////////
+//
+void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
+{
+  if(aPartDef != currentParticleDef) {
+  
+        currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
+        massRatio=1;
+        if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
+        currentParticleIndex=1000000;
+        for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
+                if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
+        }       
+         
+  }  
+}
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+//
+G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
                                 anAdjointCSMatrix,G4double Tcut)
 { 
-  std::vector< G4double > *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
+  std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
   if (theLogPrimEnergyVector->size() ==0){
-        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
-        G4cout<<"The sampling procedure will be stopped."<<std::endl;
+        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
+        G4cout<<"The s"<<G4endl;
         return 0.;
         
   }
-  //G4cout<<"A prim/Tcut "<<aPrimEnergy<<'\t'<<Tcut<<std::endl; 
   G4double log_Tcut = std::log(Tcut);
   G4double log_E =std::log(aPrimEnergy);
@@ -834 +967 @@
  
   size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
-  //G4cout<<"Prim energy "<<(*thePrimEnergyVector)[0]<<std::endl;
-  //G4cout<<"Prim energy[ind]"<<(*thePrimEnergyVector)[ind]<<std::endl;
-  //G4cout<<"Prim energy ind"<<ind<<std::endl;
-  
   G4double aLogPrimEnergy1,aLogPrimEnergy2;
   G4double aLogCS1,aLogCS2;
   G4double log01,log02;
-  std::vector< G4double>* aLogSecondEnergyVector1 =0;
-  std::vector< G4double>* aLogSecondEnergyVector2  =0;
-  std::vector< G4double>* aLogProbVector1=0;
-  std::vector< G4double>* aLogProbVector2=0;
+  std::vector< double>* aLogSecondEnergyVector1 =0;
+  std::vector< double>* aLogSecondEnergyVector2  =0;
+  std::vector< double>* aLogProbVector1=0;
+  std::vector< double>* aLogProbVector2=0;
   std::vector< size_t>* aLogProbVectorIndex1=0;
   std::vector< size_t>* aLogProbVectorIndex2=0;
@@ -851 +980 @@
   anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
   anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
-  //G4cout<<"aSecondEnergyVector1.size() "<<aSecondEnergyVector1->size()<<std::endl;
-  //G4cout<<aSecondEnergyVector1<<std::endl;
-  //G4cout<<"aSecondEnergyVector2.size() "<<aSecondEnergyVector2->size()<<std::endl;
   if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
         G4double log_minimum_prob1, log_minimum_prob2;
-        
-        //G4cout<<aSecondEnergyVector1->size()<<std::endl;
         log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
         log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
-        //G4cout<<"minimum_prob1 "<< std::exp(log_minimum_prob1)<<std::endl;
-        //G4cout<<"minimum_prob2 "<< std::exp(log_minimum_prob2)<<std::endl;
-        //G4cout<<"Tcut "<<std::endl;
         aLogCS1+= log_minimum_prob1;
         aLogCS2+= log_minimum_prob2;

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointCSMatrix.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+
 #include "G4AdjointCSMatrix.hh"
 #include <iomanip>
 #include <fstream>
-
 #include "G4AdjointInterpolator.hh"
-
 ///////////////////////////////////////////////////////
 //
@@ -64 +66 @@
 ///////////////////////////////////////////////////////
 //
- void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< G4double>* aLogSecondEnergyVector,
-                                                                       std::vector< G4double>* aLogProbVector,size_t n_pro_decade){
+ void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< double>* aLogSecondEnergyVector,
+                                                                       std::vector< double>* aLogProbVector,size_t n_pro_decade){
         
         G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
-        //Add this time we consider that the energy are given monotically
         
+        //At this time we consider that the energy is increasing monotically
         theLogPrimEnergyVector.push_back(aLogPrimEnergy);
         theLogCrossSectionVector.push_back(aLogCS);
         theLogSecondEnergyMatrix.push_back(aLogSecondEnergyVector);
-        //G4cout<<"Test Add Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
-        //G4cout<<theSecondEnergyMatrix.size()<<std::endl;
         theLogProbMatrix.push_back(aLogProbVector);
-        //G4cout<<"Test Add Data 1 "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
-        //G4cout<<theSecondEnergyMatrix.size()<<std::endl;
+        
         std::vector< size_t>* aLogProbVectorIndex = 0;
         dlog =0;
+        
         if (n_pro_decade > 0 && aLogProbVector->size()>0) {
                 aLogProbVectorIndex = new std::vector< size_t>();
@@ -85 +85 @@
                 G4double log_val = int(std::min((*aLogProbVector)[0],aLogProbVector->back())/dlog)*dlog;
                 log0Vector.push_back(log_val);
+                
                 while(log_val<0.) {
                         aLogProbVectorIndex->push_back(theInterpolator->FindPosition(log_val,(*aLogProbVector)));
@@ -102 +103 @@
 ///////////////////////////////////////////////////////
 //
-bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< G4double>*& aLogSecondEnergyVector,
-                                                                       std::vector< G4double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
+G4bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< double>*& aLogSecondEnergyVector,
+                                                                       std::vector< double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
 {       if (i>= nb_of_PrimEnergy) return false;
-        //G4cout<<"Test Get Data "<<std::endl;
+        //G4cout<<"Test Get Data "<<G4endl;
         aLogPrimEnergy = theLogPrimEnergyVector[i];
         aLogCS = theLogCrossSectionVector[i];
         aLogSecondEnergyVector = theLogSecondEnergyMatrix[i];
-        //G4cout<<"Test Get Data "<<this<<'\t'<<theSecondEnergyMatrix[i]->size()<<std::endl;
-        //G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
-        //G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector<<std::endl;
         aLogProbVector = theLogProbMatrix[i];
         aLogProbVectorIndex = theLogProbMatrixIndex[i];
         log0=log0Vector[i];
-        //G4cout<<"Test Get Data 1 "<<this<<'\t'<<theProbMatrix[i]->size()<<std::endl;
-        //G4cout<<"Test Get Data 1 "<<this<<'\t'<<aProbVector->size()<<std::endl;
-        //G4cout<<"Test Get Data 1 "<<this<<'\t'<<aLogProbVectorIndex<<std::endl;
         return true;
         
@@ -127 +122 @@
         FileOutput<<std::setiosflags(std::ios::scientific);
         FileOutput<<std::setprecision(6);
-        FileOutput<<theLogPrimEnergyVector.size()<<std::endl;
+        FileOutput<<theLogPrimEnergyVector.size()<<G4endl;
         for (size_t i=0;i<theLogPrimEnergyVector.size();i++){
-                FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<std::endl;
+                FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<G4endl;
                 size_t j1=0;
-                FileOutput<<theLogSecondEnergyMatrix[i]->size()<<std::endl;
+                FileOutput<<theLogSecondEnergyMatrix[i]->size()<<G4endl;
                 for (size_t j=0;j<theLogSecondEnergyMatrix[i]->size();j++){
                         FileOutput<<std::exp((*theLogSecondEnergyMatrix[i])[j]);
@@ -137 +132 @@
                         if (j1<10) FileOutput<<'\t';
                         else  {
-                                FileOutput<<std::endl;
+                                FileOutput<<G4endl;
                                 j1=0;
                         }       
                 }
-                if (j1>0) FileOutput<<std::endl;
+                if (j1>0) FileOutput<<G4endl;
                 j1=0;
-                FileOutput<<theLogProbMatrix[i]->size()<<std::endl;
+                FileOutput<<theLogProbMatrix[i]->size()<<G4endl;
                 for (size_t j=0;j<theLogProbMatrix[i]->size();j++){
                         FileOutput<<std::exp((*theLogProbMatrix[i])[j]);
@@ -149 +144 @@
                         if (j1<10) FileOutput<<'\t';
                         else  {
-                                FileOutput<<std::endl;
+                                FileOutput<<G4endl;
                                 j1=0;
                         }       
                 }
-                if (j1>0) FileOutput<<std::endl;
+                if (j1>0) FileOutput<<G4endl;
                 
                 
@@ -177 +172 @@
                 theLogCrossSectionVector.push_back(CS);
                 FileOutput>>n2;
-                theLogSecondEnergyMatrix.push_back(new std::vector<double>());
-                theLogProbMatrix.push_back(new std::vector<double>());
+                theLogSecondEnergyMatrix.push_back(new std::vector<G4double>());
+                theLogProbMatrix.push_back(new std::vector<G4double>());
                 
                 for (size_t j=0; j<n2;j++){

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointComptonModel.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointComptonModel.hh"
 #include "G4AdjointCSManager.hh"
@@ -34 +37 @@
 #include "G4AdjointGamma.hh"
 #include "G4Gamma.hh"
+#include "G4KleinNishinaCompton.hh"
 
 
@@ -44 +48 @@
   SetUseMatrix(true);
   SetUseMatrixPerElement(true);
-  SetIsIonisation(false);
   SetUseOnlyOneMatrixForAllElements(true);
   theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
@@ -50 +53 @@
   theDirectPrimaryPartDef=G4Gamma::Gamma();
   second_part_of_same_type=false;
+  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
  
 }
@@ -62 +66 @@
                        G4ParticleChange* fParticleChange)
 { 
-   
+   if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
    
    //A recall of the compton scattering law is 
@@ -71 +75 @@
    
   const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
-  //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
-  size_t ind= 0;
-  
-  
- 
- 
- //Elastic inverse scattering //not correct in all the cases 
- //---------------------------------------------------------
- G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
- 
- //G4cout<<adjointPrimKinEnergy<<std::endl;
- if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
+  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
+  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
         return;
- }
+  }
+ 
  
  //Sample secondary energy
  //-----------------------
  G4double gammaE1;
-
- gammaE1 = SampleAdjSecEnergyFromCSMatrix(ind,
-                                                  adjointPrimKinEnergy,
+ gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
                                                   IsScatProjToProjCase);
  
@@ -108 +101 @@
  //Cos th
  //-------
-// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<std::endl;
+// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl;
  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
  if (!IsScatProjToProjCase) {
@@ -116 +109 @@
  G4double sin_th = 0.;
  if (std::abs(cos_th)>1){
-        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<std::endl;
+        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
         if (cos_th>0) {
                 cos_th=1.;
@@ -136 +129 @@
  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
  gammaMomentum1.rotateUz(dir_parallel);
-// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<std::endl;
+// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl;
  
  
  //It is important to correct the weight of particles before adding the secondary
  //------------------------------------------------------------------------------
- CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,gammaE1);
- 
- if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
+ CorrectPostStepWeight(fParticleChange, 
+                        aTrack.GetWeight(), 
+                        adjointPrimKinEnergy,
+                        gammaE1,
+                        IsScatProjToProjCase);
+ 
+ if (!IsScatProjToProjCase){ //kill the primary and add a secondary
         fParticleChange->ProposeTrackStatus(fStopAndKill);
         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
-        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<std::endl;
+        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
  }
  else {
@@ -154 +151 @@
  
         
-}                       
+} 
+////////////////////////////////////////////////////////////////////////////////
+//
+void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
+                       G4bool IsScatProjToProjCase,
+                       G4ParticleChange* fParticleChange)
+{ 
+
+ const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
+ DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
+ 
+ 
+ G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
+
+ 
+ if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
+        return;
+ }
+  
+ 
+ G4double diffCSUsed=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2; 
+ G4double gammaE1=0.;
+ G4double gammaE2=0.;
+ if (!IsScatProjToProjCase){
+        
+        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
+        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
+        if (Emin>=Emax) return;
+        G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
+        G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
+        gammaE1=adjointPrimKinEnergy/(1.-f1*pow(f2,G4UniformRand()));;
+        gammaE2=gammaE1-adjointPrimKinEnergy;
+        diffCSUsed= diffCSUsed*(1.+2.*log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
+        
+        
+ }
+ else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
+        G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
+        if (Emin>=Emax) return;
+        gammaE2 =adjointPrimKinEnergy;
+        gammaE1=Emin*pow(Emax/Emin,G4UniformRand());
+        diffCSUsed= diffCSUsed/gammaE1;
+        
+ }
+  
+  
+  
+                                                        
+ //Weight correction
+ //-----------------------
+ //First w_corr is set to the ratio between adjoint total CS and fwd total CS
+ G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
+
+ //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
+ 
+ 
+ G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
+ if (diffCS >0)  diffCS /=G4direct_CS;  // here we have the normalised diffCS
+ diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
+ //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
+ //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;                                 
+ 
+ w_corr*=diffCS/diffCSUsed;
+           
+ G4double new_weight = aTrack.GetWeight()*w_corr;
+ fParticleChange->SetParentWeightByProcess(false);
+ fParticleChange->SetSecondaryWeightByProcess(false);
+ fParticleChange->ProposeParentWeight(new_weight); 
+ 
+ 
+ 
+ //Cos th
+ //-------
+
+ G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
+ if (!IsScatProjToProjCase) {
+        G4double p_elec=theAdjointPrimary->GetTotalMomentum();
+        cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
+ }
+ G4double sin_th = 0.;
+ if (std::abs(cos_th)>1){
+        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
+        if (cos_th>0) {
+                cos_th=1.;
+        }
+        else    cos_th=-1.;
+        sin_th=0.;
+ }
+ else  sin_th = std::sqrt(1.-cos_th*cos_th);
+
+ 
+ 
+ 
+ //gamma0 momentum
+ //--------------------
+
+ 
+ G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
+ G4double phi =G4UniformRand()*2.*3.1415926;
+ G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
+ gammaMomentum1.rotateUz(dir_parallel);
+ 
+ 
+
+ 
+ if (!IsScatProjToProjCase){ //kill the primary and add a secondary
+        fParticleChange->ProposeTrackStatus(fStopAndKill);
+        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
+        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
+ }
+ else {
+        fParticleChange->ProposeEnergy(gammaE1);
+        fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
+ }
+ 
+ 
+ 
+} 
+
+                        
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -179 +296 @@
                                       G4double )
 { //Based on Klein Nishina formula
- //* In the forward case (see G4KleinNishinaModel)  the  cross section is parametrised while the secondaries are sampled from the 
+ // In the forward case (see G4KleinNishinaModel)  the  cross section is parametrised while
+ // the secondaries are sampled from the 
  // Klein Nishida differential cross section
- // The used diffrential cross section here is therefore the cross section multiplied by the normalidsed differential Klein Nishida cross section
+ // The used diffrential cross section here is therefore the cross section multiplied by the normalised 
+ //differential Klein Nishida cross section
  
  
@@ -191 +310 @@
  G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
  if (gamEnergy1 >gamEnergy1_max ||  gamEnergy1<gamEnergy1_min) {
-        /*G4cout<<"the differential CS is null"<<std::endl;
-        G4cout<<gamEnergy0<<std::endl;
-        G4cout<<gamEnergy1<<std::endl;
-        G4cout<<gamEnergy1_min<<std::endl;*/
+        /*G4cout<<"the differential CS is null"<<G4endl;
+        G4cout<<gamEnergy0<<G4endl;
+        G4cout<<gamEnergy1<<G4endl;
+        G4cout<<gamEnergy1_min<<G4endl;*/
         return 0.;
  }
@@ -222 +341 @@
  //-------------------------------
  
- G4double G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
+ G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
                                              gamEnergy0,
                                              Z, 0., 0.,0.);
  
  dCS_dE1 *= G4direct_CS/CS;
-/* G4cout<<"the differential CS is not null"<<std::endl;
- G4cout<<gamEnergy0<<std::endl;
- G4cout<<gamEnergy1<<std::endl;*/
+/* G4cout<<"the differential CS is not null"<<G4endl;
+ G4cout<<gamEnergy0<<G4endl;
+ G4cout<<gamEnergy1<<G4endl;*/
  
  return dCS_dE1;
@@ -235 +354 @@
 
 }
+                                      
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -251 +371 @@
   return  emin; 
 }
+////////////////////////////////////////////////////////////////////////////////
+//
+G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
+                                             G4double primEnergy,
+                                             G4bool IsScatProjToProjCase)
+{ 
+  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
+  DefineCurrentMaterial(aCouple);
+  
+  
+  G4double Cross=0.;
+  G4double Emax_proj =0.;
+  G4double Emin_proj =0.;
+  if (!IsScatProjToProjCase ){
+        Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
+        Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
+        if (Emax_proj>Emin_proj ){
+                 Cross= log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
+                                                                *(1.+2.*log(1.+electron_mass_c2/primEnergy));
+        }        
+  }
+  else {
+        Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
+        Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
+        if (Emax_proj>Emin_proj) {
+                Cross = log(Emax_proj/Emin_proj);
+                //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
+        }
+        
+        
+  }
+  
+  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
+  lastCS=Cross;
+  return Cross; 
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointInterpolator.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointCSMatrix.hh"
 #include "G4AdjointInterpolator.hh"
@@ -61 +64 @@
 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
 { G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
-  //G4cout<<"Linear "<<res<<std::endl;
+  //G4cout<<"Linear "<<res<<G4endl;
   return res;   
 }
@@ -69 +72 @@
 { if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
   G4double B=std::log(y2/y1)/std::log(x2/x1);
-  //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<std::endl;
+  //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<G4endl;
   G4double A=y1/std::pow(x1,B);
   G4double res=A*std::pow(x,B);
- // G4cout<<"Log "<<res<<std::endl;
+ // G4cout<<"Log "<<res<<G4endl;
   return res;
 }
@@ -98 +101 @@
   }
   else {
-        //G4cout<<"The interpolation method that you invoked does not exist!"<<std::endl; 
+        //G4cout<<"The interpolation method that you invoked does not exist!"<<G4endl; 
         return -1111111111.;
   }
@@ -104 +107 @@
 ///////////////////////////////////////////////////////
 //
-size_t  G4AdjointInterpolator::FindPosition(G4double& x,std::vector<double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
+size_t  G4AdjointInterpolator::FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
 {  //most rapid nethod could be used probably
-   //It is important to put std::vector<double>& such that the vector itself is used and not a copy
+   //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
   
   
@@ -149 +152 @@
 ///////////////////////////////////////////////////////
 //
-size_t  G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<double>& log_x_vec) //only valid if x_vec is monotically increasing
+size_t  G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec) //only valid if x_vec is monotically increasing
 {  //most rapid nethod could be used probably
-   //It is important to put std::vector<double>& such that the vector itself is used and not a copy
-  
+   //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
+  return FindPosition(log_x, log_x_vec);
   if (log_x_vec.size()>3){ 
         size_t ind=0;
@@ -170 +173 @@
 ///////////////////////////////////////////////////////
 //
-G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,G4String InterPolMethod)
+G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
 { size_t i=FindPosition(x,x_vec);
-  //G4cout<<i<<std::endl;
-  //G4cout<<x<<std::endl;
-  //G4cout<<x_vec[i]<<std::endl; 
+  //G4cout<<i<<G4endl;
+  //G4cout<<x<<G4endl;
+  //G4cout<<x_vec[i]<<G4endl; 
   return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
 }
@@ -180 +183 @@
 ///////////////////////////////////////////////////////
 //
-G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,
+G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
                                             std::vector<size_t>& index_vec,G4double x0, G4double dx) //only linear interpolation possible
 { size_t ind=0;
@@ -202 +205 @@
 ///////////////////////////////////////////////////////
 //
-G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<double>& log_x_vec,std::vector<double>& log_y_vec)
+G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
 { //size_t i=0;
   size_t i=FindPositionForLogVector(log_x,log_x_vec);
-  /*G4cout<<"In interpolate "<<std::endl;
-  G4cout<<i<<std::endl;
-  G4cout<<log_x<<std::endl;
-  G4cout<<log_x_vec[i]<<std::endl;
-  G4cout<<log_x_vec[i+1]<<std::endl;
-  G4cout<<log_y_vec[i]<<std::endl;
-  G4cout<<log_y_vec[i+1]<<std::endl;*/
+  /*G4cout<<"In interpolate "<<G4endl;
+  G4cout<<i<<G4endl;
+  G4cout<<log_x<<G4endl;
+  G4cout<<log_x_vec[i]<<G4endl;
+  G4cout<<log_x_vec[i+1]<<G4endl;
+  G4cout<<log_y_vec[i]<<G4endl;
+  G4cout<<log_y_vec[i+1]<<G4endl;*/
   
   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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4AdjointPhotoElectricModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4AdjointPhotoElectricModel.hh"
 #include "G4AdjointCSManager.hh"
@@ -35 +38 @@
 #include "G4AdjointGamma.hh"
 
+
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -41 +45 @@
 
 { SetUseMatrix(false);
+  SetApplyCutInRange(false);
   current_eEnergy =0.;
   totAdjointCS=0.;
+  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
+  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
+  theDirectPrimaryPartDef=G4Gamma::Gamma();
+  second_part_of_same_type=false;
+  theDirectPEEffectModel = new G4PEEffectModel();
+ 
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -57 +68 @@
 
   //Compute the totAdjointCS vectors if not already done for the current couple and electron energy
+  //-----------------------------------------------------------------------------------------------
   const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
   const G4DynamicParticle* aDynPart =  aTrack.GetDynamicParticle() ;
   G4double electronEnergy = aDynPart->GetKineticEnergy();
   G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
-  totAdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
+  pre_step_AdjointCS = totAdjointCS; //The last computed CS was  at pre step point
+  G4double adjCS;
+  adjCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
+  post_step_AdjointCS = totAdjointCS; 
                                 
- 
-  //Sample gamma energy
-  //-------------
-        /////////////////////////////////////////////////////////////////////////////////
-//      Module:         G4ContinuousGainOfEnergy.hh
-//      Author:         L. Desorgher
-//      Date:           1 September 2007
-//      Organisation:   SpaceIT GmbH
-//      Customer:       ESA/ESTEC
-/////////////////////////////////////////////////////////////////////////////////
-//
-// CHANGE HISTORY
-// --------------
-//      ChangeHistory: 
-//              1 September 2007 creation by L. Desorgher               
-//
-//-------------------------------------------------------------
-//      Documentation:
-//              Modell for the adjoint compton scattering
-//
+
+
 
   //Sample element
   //-------------
    const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
-   const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
    size_t nelm =  currentMaterial->GetNumberOfElements();
-   G4double rand_CS= totAdjointCS*G4UniformRand();
+   G4double rand_CS= G4UniformRand()*xsec[nelm-1];
    for (index_element=0; index_element<nelm-1; index_element++){
         if (rand_CS<xsec[index_element]) break;
@@ -96 +92 @@
    //Sample shell and binding energy
    //-------------
-   rand_CS= totAdjointCS*G4UniformRand()/theAtomNumDensityVector[index_element];
    G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
+   rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand();
    G4int i  = 0;  
    for (i=0; i<nShells-1; i++){
@@ -113 +109 @@
    G4double gamma   = 1. + electronEnergy/electron_mass_c2;
    if (gamma <= 5.) {
-        G4double beta  = std::sqrt(gamma*gamma-1.)/gamma;
+        G4double beta  = sqrt(gamma*gamma-1.)/gamma;
         G4double b     = 0.5*gamma*(gamma-1.)*(gamma-2);
 
@@ -131 +127 @@
   
      
-  G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
+  G4double sin_theta = sqrt(1.-cos_theta*cos_theta);
   G4double Phi     = twopi * G4UniformRand();
-  G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
+  G4double dirx = sin_theta*cos(Phi),diry = sin_theta*sin(Phi),dirz = cos_theta;
   G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz);
   adjoint_gammaDirection.rotateUz(electronDirection);
@@ -141 +137 @@
   //Weight correction
  //-----------------------                                         
-  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy);       
+  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase);  
  
   
@@ -149 +145 @@
   G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
                        G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
+  
+  
+  
+ 
+  
   fParticleChange->ProposeTrackStatus(fStopAndKill);
-  fParticleChange->AddSecondary(anAdjointGamma);
-     
+  fParticleChange->AddSecondary(anAdjointGamma);    
         
 
   
 
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+void G4AdjointPhotoElectricModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, 
+                                                            G4double old_weight,  
+                                                            G4double adjointPrimKinEnergy, 
+                                                            G4double projectileKinEnergy ,
+                                                            G4bool  ) 
+{
+ G4double new_weight=old_weight;
+
+ G4double w_corr =G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection()/factorCSBiasing;
+ w_corr*=post_step_AdjointCS/pre_step_AdjointCS; 
+ new_weight*=w_corr; 
+ new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
+ fParticleChange->SetParentWeightByProcess(false);
+ fParticleChange->SetSecondaryWeightByProcess(false);
+ fParticleChange->ProposeParentWeight(new_weight);      
 }       
 
@@ -164 +183 @@
                                 G4double electronEnergy,
                                 G4bool IsScatProjToProjCase)
-{ if (IsScatProjToProjCase) return 0.;
+{ 
+  
+
+  if (IsScatProjToProjCase)  return 0.;
+
+        
   if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
         totAdjointCS = 0.;
         DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
         const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
-        const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
+        const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
         size_t nelm =  currentMaterial->GetNumberOfElements();
         for (index_element=0;index_element<nelm;index_element++){
@@ -176 +200 @@
                 xsec[index_element] = totAdjointCS;
         } 
-  }     
-  return totAdjointCS;
-  
+
+        totBiasedAdjointCS=std::min(totAdjointCS,0.01);
+//      totBiasedAdjointCS=totAdjointCS;
+        factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
+        lastCS=totBiasedAdjointCS;
+        
+        
+  }
+  return totBiasedAdjointCS;
+
   
 }
@@ -188 +219 @@
   G4int nShells = anElement->GetNbOfAtomicShells();
   G4double Z= anElement->GetZ();
-  G4double N= anElement->GetN();
   G4int i  = 0;  
   G4double B0=anElement->GetAtomicShell(0);
   G4double gammaEnergy = electronEnergy+B0;
-  G4double adjointCS = theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,N,0.,0.)*electronEnergy/gammaEnergy; 
+  G4double CS= theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
+  G4double adjointCS =0.;
+  if (CS >0) adjointCS += CS/gammaEnergy; 
   shell_prob[index_element][0] = adjointCS;                                          
   for (i=1;i<nShells;i++){
-        //G4cout<<i<<std::endl;
+        //G4cout<<i<<G4endl;
         G4double Bi_= anElement->GetAtomicShell(i-1);
         G4double Bi = anElement->GetAtomicShell(i);
-        //G4cout<<Bi_<<'\t'<<Bi<<std::endl;
+        //G4cout<<Bi_<<'\t'<<Bi<<G4endl;
         if (electronEnergy <Bi_-Bi) {
                 gammaEnergy = electronEnergy+Bi;
-                adjointCS +=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,anElement->GetZ(),N,0.,0.)*electronEnergy/gammaEnergy;
+                
+                CS=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
+                if (CS>0) adjointCS +=CS/gammaEnergy;
         }
         shell_prob[index_element][i] = adjointCS;       
   
   }
-  
+  adjointCS*=electronEnergy;
   return adjointCS;
   

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4ContinuousGainOfEnergy.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4ContinuousGainOfEnergy.hh"
 #include "G4Step.hh"
@@ -31 +34 @@
 #include "G4VParticleChange.hh"
 #include "G4UnitsTable.hh"
+#include "G4AdjointCSManager.hh"
+#include "G4LossTableManager.hh"
 
 
@@ -38 +43 @@
   G4ProcessType type): G4VContinuousProcess(name, type)
 {
+
 
   linLossLimit=0.05;
@@ -44 +50 @@
   is_integral = false;
   
+  //Will be properly set in SetDirectParticle()
+  IsIon=false;
+  massRatio =1.;
+  chargeSqRatio=1.;
+  preStepChargeSqRatio=1.;
+
+  
+
+  
+  
 }
 
@@ -70 +86 @@
 }
 
-
-
+///////////////////////////////////////////////////////
+//
+void  G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
+{theDirectPartDef=p;
+ if (theDirectPartDef->GetParticleType()== "nucleus") {
+         IsIon=true;
+         massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass();
+         G4double q=theDirectPartDef->GetPDGCharge();
+         chargeSqRatio=q*q;
+        
+         
+ }
+ 
+}
 
 ///////////////////////////////////////////////////////
@@ -80 +108 @@
 {
    
+  //Caution in this method the  step length should be the true step length
+  // 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 
+  //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
+  //
+  
+  
+  
   aParticleChange.Initialize(track);
   
@@ -86 +121 @@
   G4double length = step.GetStepLength();
   G4double degain  = 0.0;
- 
- 
- 
-
+  
+  
+ 
   // Compute this for weight change after continuous energy loss
   //-------------------------------------------------------------
-  G4double DEDX_before = 
-                theDirectEnergyLossProcess
-                                ->GetDEDX(preStepKinEnergy, currentCouple);
- 
- 
+  G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple);
+   
+  
   
   // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
@@ -103 +135 @@
   G4DynamicParticle* dynParticle = new G4DynamicParticle();
   *dynParticle = *(track.GetDynamicParticle());
-  G4double Tkin = dynParticle->GetKineticEnergy();
-  
+  dynParticle->SetDefinition(theDirectPartDef);
+  G4double Tkin = dynParticle->GetKineticEnergy(); 
+  G4double Tkin1=Tkin*0.001;
+
   size_t n=1;
   if (is_integral ) n=10;
+  n=1;
   G4double dlength= length/n; 
   for (size_t i=0;i<n;i++) {
+        G4double factor_dE=1.;
+        if (Tkin != preStepKinEnergy && IsIon) {
+                chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
+                theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio); 
+        
+        }
+  
         G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple);
         if( dlength <= linLossLimit * r ) {
                 degain = DEDX_before*dlength;
+                G4double degain1 = dlength*theDirectEnergyLossProcess->GetDEDX(Tkin1, currentCouple);
+                factor_dE=1.+(degain1-degain)/(Tkin1-Tkin);
         } 
         else {
-                G4double x = r + length;
-                degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
+                G4double x = r + dlength;
+                //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
+                G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
+                if (IsIon){
+                        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
+                        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
+                        G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
+                        while (std::abs(x-x1)>0.01*x) {
+                                E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
+                                chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
+                                theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
+                                x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
+                        
+                        } 
+                }
+                G4double r1 = theDirectEnergyLossProcess->GetRange(Tkin1, currentCouple);
+                G4double x1 = r1 + dlength;
+                G4double E1 = theDirectEnergyLossProcess->GetKineticEnergy(x1,currentCouple);
+                factor_dE=(E1-E)/(Tkin1-Tkin);
+                degain=E-Tkin;  
+                
+                
+                
         }
-        G4VEmModel* currentModel =  theDirectEnergyLossProcess->SelectModelForMaterial(Tkin+degain,currentMaterialIndex);
+        //G4cout<<degain<<G4endl;
         G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
         tmax = std::min(tmax,currentTcut);
- 
-  
+        
+        
+        dynParticle->SetKineticEnergy(Tkin+degain);
+
+        // Corrections, which cannot be tabulated for ions
+        //----------------------------------------
+        G4double esecdep=0;//not used in most models
+        currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength); 
 
         // Sample fluctuations
         //-------------------
+        
         
         G4double deltaE =0.;
         if (lossFluctuationFlag ) {
                 deltaE = currentModel->GetModelOfFluctuations()->
-                                                SampleFluctuations(currentMaterial,dynParticle,tmax,length,degain)-degain;
+                                                SampleFluctuations(currentMaterial,dynParticle,tmax,dlength,degain)-degain;
         }
-        Tkin+=degain+deltaE;
+        
+        G4double egain=degain+deltaE;
+        if (egain <=0) egain=degain;
+        Tkin+=egain;
         dynParticle->SetKineticEnergy(Tkin);
-        
+ }
+ 
+ 
+  
+
+  
+  delete dynParticle;
+ 
+  if (IsIon){
+        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
+        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
+                
   }
- 
-  
-
-  // Corrections, which cannot be tabulated 
-  // probably  this should be also changed
-  // at this time it does nothing so we can leave it
-  //CorrectionsAlongStep(currentCouple, dynParticle, egain, length);
-  
-  delete dynParticle;
- 
   
   G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
-  G4double weight_correction=DEDX_after/DEDX_before; //probably not needed
-  weight_correction=1.;
+  
+  
+  G4double weight_correction=DEDX_after/DEDX_before;
+ 
   
   aParticleChange.ProposeEnergy(Tkin);
@@ -154 +232 @@
   //we still need to register in the particleChange the modification of the weight of the particle 
   G4double new_weight=weight_correction*track.GetWeight();
-  aParticleChange.SetParentWeightByProcess(true);
+  aParticleChange.SetParentWeightByProcess(false);
   aParticleChange.ProposeParentWeight(new_weight);
   
@@ -168 +246 @@
   lossFluctuationFlag = val;
 }
+///////////////////////////////////////////////////////
+//
+
+
+
+G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
+                G4double , G4double , G4double& )
+{ 
+  G4double x = DBL_MAX;
+  x=.1*mm;
+ 
+ 
+  DefineMaterial(track.GetMaterialCutsCouple());
+ 
+  preStepKinEnergy = track.GetKineticEnergy(); 
+  preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
+  currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
+  G4double emax_model=currentModel->HighEnergyLimit();
+  if (IsIon) {
+        chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
+        preStepChargeSqRatio = chargeSqRatio;
+        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
+  } 
+  
+  
+  G4double maxE =1.1*preStepKinEnergy;
+  /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
+  else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
+  else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
+   
+  if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
+ 
+  maxE=std::min(emax_model*1.001,maxE);
+        
+  G4double r = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
+  
+  if (IsIon) {
+        G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
+        theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
+  }     
+  
+  G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
+  
+  if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
+  
+  
+
+  x=r1-r;
+  x=std::max(r1-r,0.001*mm);
+ 
+  return x;
+  
+ 
+}
+#include "G4EmCorrections.hh"
+///////////////////////////////////////////////////////
+//
+
+void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
+{ 
+
+  G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy); 
+  if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4InversePEEffect.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4InversePEEffect.hh"
 #include "G4VEmAdjointModel.hh"
@@ -31 +34 @@
 //
 G4InversePEEffect::G4InversePEEffect(G4String process_name,G4AdjointPhotoElectricModel* aModel):
-                                G4VAdjointInverseScattering(process_name,false)
+                                G4VAdjointReverseReaction(process_name,false)
 {theAdjointEMModel = aModel;
- theAdjointEMModel->SetSecondPartOfSameType(false); 
+ theAdjointEMModel->SetSecondPartOfSameType(false);
+ SetIntegralMode(false);
+  
+ 
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4VEmAdjointModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4VEmAdjointModel.hh"
 #include "G4AdjointCSManager.hh"
-
-
 #include "G4Integrator.hh"
 #include "G4TrackStatus.hh"
@@ -33 +34 @@
 #include "G4AdjointElectron.hh"
 #include "G4AdjointInterpolator.hh"
+#include "G4PhysicsTable.hh"
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -39 +41 @@
 name(nam)
 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
-{ G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
-  CorrectWeightMode =true;
-  UseMatrix =true;
-  UseMatrixPerElement = true;
-  ApplyCutInRange = true;
-  ApplyBiasing = true;
-  UseOnlyOneMatrixForAllElements = true;
-  IsIonisation =true; 
-  CS_biasing_factor =1.;
-  //ApplyBiasing = false;
+{ 
+  G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
+  second_part_of_same_type =false;
+  theDirectEMModel=0;
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -54 +50 @@
 G4VEmAdjointModel::~G4VEmAdjointModel()
 {;}
-////////////////////////////////////////////////////////////////////////////////
-//
-void G4VEmAdjointModel::SampleSecondaries(const G4Track& aTrack,
-                       G4bool IsScatProjToProjCase,
-                       G4ParticleChange* fParticleChange)
-{ 
-
-  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
-  //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
-  size_t ind=0;
-  if (!UseMatrixPerElement) ind = currentMaterialIndex;
-  //G4cout<<theAdjointPrimary<<std::endl;
-  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
-        std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
-        if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
-        G4double rand_var= G4UniformRand();
-        G4double SumCS=0.;
-        for (size_t i=0;i<CS_Vs_Element->size();i++){
-                SumCS+=(*CS_Vs_Element)[i];
-                if (rand_var<=SumCS/lastCS){
-                        ind=i;
-                        break;
-                }
-        }
-        ind = currentMaterial->GetElement(ind)->GetIndex();
-  }
-  
- 
- 
- //Elastic inverse scattering //not correct in all the cases 
- //---------------------------------------------------------
- G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
- G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
- G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
- //G4cout<<adjointPrimKinEnergy<<std::endl;
- if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
-        return;
- }
- 
- //Sample secondary energy
- //-----------------------
- G4double projectileKinEnergy;
-// if (!IsIonisation  ) {
-        projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
-                                                  adjointPrimKinEnergy,
-                                                  IsScatProjToProjCase);
- //}
- /*else {
-        projectileKinEnergy =   SampleAdjSecEnergyFromDiffCrossSectionPerAtom(adjointPrimKinEnergy,IsScatProjToProjCase);
-        //G4cout<<projectileKinEnergy<<std::endl;
- }*/    
- //Weight correction
- //-----------------------
- 
- CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);                                     
- 
-                 
- 
- 
- 
- //Kinematic
- //---------
- 
- G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
- G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
- G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
- 
- 
- 
- //Companion
- //-----------
- G4double companionM0;
- companionM0=(adjointPrimTotalEnergy-adjointPrimKinEnergy);
- if (IsScatProjToProjCase) {
-        companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
- } 
- G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
- G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
- 
- 
- //Projectile momentum
- //--------------------
- G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
- G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
- G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
- G4double phi =G4UniformRand()*2.*3.1415926;
- G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
- projectileMomentum.rotateUz(dir_parallel);
- 
- 
- 
- if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
-        fParticleChange->ProposeTrackStatus(fStopAndKill);
-        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
-        //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
- }
- else {
-        fParticleChange->ProposeEnergy(projectileKinEnergy);
-        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
- }      
-}
-////////////////////////////////////////////////////////////////////////////////
-//
-void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,  G4double , G4double ) 
-{
- G4double new_weight=old_weight;
- if (CorrectWeightMode) {
-        G4double w_corr =1./CS_biasing_factor;
-        //G4cout<<w_corr<<std::endl;
-        
-        /*G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(theAdjEquivOfDirectPrimPartDef, 
-                                                                               theAdjEquivOfDirectSecondPartDef,
-                                                                               adjointPrimKinEnergy,projectileKinEnergy,
-                                                                               aTrack.GetMaterialCutsCouple());
-        w_corr = projectileKinEnergy;
-        G4double Emin,Emax;
-        if (IsScatProjToProjCase) {
-                Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
-                Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy, currentTcutForDirectSecond);
-                
-        }
-        else {
-                Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
-                Emin = GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);    
-        }
-        w_corr *=std::log(Emax/Emin)/(Emax-Emin); */    
-                
-        new_weight*=w_corr;
- }
- G4cout<< "new weight"<<new_weight<<std::endl;
- fParticleChange->SetParentWeightByProcess(false);
- fParticleChange->SetSecondaryWeightByProcess(false);
- fParticleChange->ProposeParentWeight(new_weight);      
-}
 ////////////////////////////////////////////////////////////////////////////////
 //                              
@@ -195 +57 @@
 { 
   DefineCurrentMaterial(aCouple);
-  //G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(G4AdjointElectron::AdjointElectron(),primEnergy,aCouple);
-  //G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(G4AdjointElectron::AdjointElectron(), primEnergy,aCouple);
-  if (IsScatProjToProjCase){
-        lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
+  preStepEnergy=primEnergy;
+  
+  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
+  if (IsScatProjToProjCase)   CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
+  lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
                                                                         this, 
                                                                         primEnergy,
                                                                         currentTcutForDirectSecond,
-                                                                        true,
-                                                                        CS_Vs_ElementForScatProjToProjCase);
-        /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectPrimPartDef,primEnergy,aCouple);
-        G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectPrimPartDef, primEnergy,aCouple);
-        */
-        //if (adjCS >0 )lastCS *=fwdCS/adjCS;
- 
-  }
-  else {
-        lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
-                                                                        this, 
-                                                                        primEnergy,
-                                                                        currentTcutForDirectSecond,
-                                                                        false,
-                                                                        CS_Vs_ElementForProdToProjCase);
-        /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectSecondPartDef,primEnergy,aCouple);
-        G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectSecondPartDef, primEnergy,aCouple);
-        */
-        //if (adjCS >0 )lastCS *=fwdCS/adjCS;                                                           
-        //lastCS=0.;                                                            
-  }
- 
-  
+                                                                        IsScatProjToProjCase,
+                                                                        *CS_Vs_Element);
+  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
+  else lastAdjointCSForProdToProjCase =lastCS;                                                          
+        
+ 
+ 
   return lastCS;
                                                                         
@@ -230 +77 @@
 ////////////////////////////////////////////////////////////////////////////////
 //
-//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine  
+//General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine  
 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
                                       G4double kinEnergyProj, 
@@ -245 +92 @@
         G4double Tmax=kinEnergyProj;
         if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
-        return Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
-        //it could be thta Tmax here should be DBLMAX
-        //Tmax=DBLMAX;
-        
-        G4double E1=kinEnergyProd;
+        
+        G4double E1=kinEnergyProd;
         G4double E2=kinEnergyProd*1.000001;
         G4double dE=(E2-E1);
@@ -256 +100 @@
         
         dSigmadEprod=(sigma1-sigma2)/dE;
-        if (dSigmadEprod>1.) {
-                G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<std::endl;
-                G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<std::endl;
-                G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<std::endl;
-                
-        }
-        
-        
-        
  }
- 
- 
-
-        
  return dSigmadEprod;   
  
@@ -307 +138 @@
         G4double Tmax=kinEnergyProj;
         if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
-        //it could be thta Tmax here should be DBLMAX
-        //Tmax=DBLMAX;
-        
-        G4double E1=kinEnergyProd;
-        
-        G4double E2=kinEnergyProd*1.0001;
+        G4double E1=kinEnergyProd;
+        G4double E2=kinEnergyProd*1.0001;
         G4double dE=(E2-E1);
         G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
-        
-        //G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
-        dSigmadEprod=sigma1/dE;
-        if (dSigmadEprod <0) { //could happen with bremstrahlung dur to suppression effect
-                G4cout<<"Halllllllllllllllllllllllllllllllllllllllllllllllo "<<kinEnergyProj<<'\t'<<E1<<'\t'<<dSigmadEprod<<std::endl;
-                E1=kinEnergyProd;
-                E2=E1*1.1;
-                dE=E2-E1;
-                sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e50);
-                G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
-                dSigmadEprod=(sigma1-sigma2)/dE;
-                G4cout<<dSigmadEprod<<std::endl;
-        }       
-        
-        
+        G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
+        dSigmadEprod=(sigma1-sigma2)/dE;
  }
-
-
-        
  return dSigmadEprod;   
  
@@ -356 +167 @@
 //
 G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
-  //return kinEnergyProj*kinEnergyProj;
-  //ApplyBiasing=false;
+  
+  
   G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;   
-  if (!ApplyBiasing) bias_factor =CS_biasing_factor; 
-  //G4cout<<bias_factor<<std::endl;
+
+
   if (UseMatrixPerElement ) {
         return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
   }
-  else {
+  else  {
         return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
   }     
 }
-//////////////////////////////////////////////////////////////////////////////
-//
-G4double G4VEmAdjointModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
-  G4double electron_mass_c2=0.51099906*MeV;
-  G4double energy = kinEnergyProj + electron_mass_c2;
-  G4double x   = kinEnergyProd/kinEnergyProj;
-  G4double gam    = energy/electron_mass_c2;
-  G4double gamma2 = gam*gam;
-  G4double beta2  = 1.0 - 1.0/gamma2;
-  
-  G4double g = (2.0*gam - 1.0)/gamma2;
-  G4double y = 1.0 - x;
-  G4double fac=twopi_mc2_rcl2/electron_mass_c2;
-  G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
-  return dCS/kinEnergyProj;
-  
- 
-
-}  
+
 ////////////////////////////////////////////////////////////////////////////////
 //
 G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
-  //return kinEnergyProj*kinEnergyProj; 
-  G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
-  //ApplyBiasing=false;
- if (!ApplyBiasing) bias_factor = CS_biasing_factor;
- //G4cout<<bias_factor<<std::endl; 
+  
+ G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
  if (UseMatrixPerElement ) {
         return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
  }      
- else {
+ else {  
         return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
  
  }      
                 
+}
+////////////////////////////////////////////////////////////////////////////////
+//
+
+G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
+{
+  return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -411 +208 @@
                                 G4double A ,
                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
-{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
-  ASelectedNucleus= G4int(A);
-  ZSelectedNucleus=G4int(Z);
+{ 
+  G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
+  ASelectedNucleus= int(A);
+  ZSelectedNucleus=int(Z);
   kinEnergyProdForIntegration = kinEnergyProd;
   
@@ -422 +220 @@
   G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
   G4double E1=minEProj;
-  std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
-  std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
+  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
+  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
   log_ESec_vector->clear();
   log_Prob_vector->clear();
@@ -429 +227 @@
   log_Prob_vector->push_back(-50.);
   
-  G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
+  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
   G4double fE=std::pow(10.,1./nbin_pro_decade);
   G4double int_cross_section=0.;
@@ -436 +234 @@
   
   while (E1 <maxEProj*0.9999999){
-        //G4cout<<E1<<'\t'<<E2<<std::endl;
-        
-        int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
-        //G4cout<<"int_cross_section 1 "<<'\t'<<int_cross_section<<std::endl;
+        //G4cout<<E1<<'\t'<<E2<<G4endl;
+        
+        int_cross_section +=integral.Simpson(this,
+        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
         log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
         log_Prob_vector->push_back(std::log(int_cross_section));        
@@ -463 +261 @@
                                 G4double A ,
                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
-{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
-  ASelectedNucleus=G4int(A);
-  ZSelectedNucleus=G4int(Z);
+{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
+  ASelectedNucleus=int(A);
+  ZSelectedNucleus=int(Z);
   kinEnergyScatProjForIntegration = kinEnergyScatProj;
   
@@ -479 +277 @@
   
   
-  std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
-  std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
+  std::vector< double>*  log_ESec_vector = new std::vector< double>();
+  std::vector< double>*  log_Prob_vector = new std::vector< double>();
   log_ESec_vector->push_back(std::log(dEmin));
   log_Prob_vector->push_back(-50.);
-  G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
+  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
   G4double fE=std::pow(dEmax/dEmin,1./nbins);
+  
+  
+  
+  
   
   G4double int_cross_section=0.;
@@ -491 +293 @@
         dE2=dE1*fE;
         int_cross_section +=integral.Simpson(this,
-        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
-        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
-        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
+        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
+        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
+        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
         log_Prob_vector->push_back(std::log(int_cross_section));        
         dE1=dE2;
 
   }
-  /*G4cout<<"total int_cross_section"<<'\t'<<int_cross_section<<std::endl;
-  G4cout<<"energy "<<kinEnergyScatProj<<std::endl;*/
-  
-  
   
   
@@ -519 +317 @@
                                 G4double kinEnergyProd,
                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
-{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
+{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
   SelectedMaterial= aMaterial;
   kinEnergyProdForIntegration = kinEnergyProd;
-  //G4cout<<aMaterial->GetName()<<std::endl;
-  //G4cout<<kinEnergyProd/MeV<<std::endl;
-  //compute the vector of integrated cross sections
+   //compute the vector of integrated cross sections
   //-------------------
   
@@ -530 +326 @@
   G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
   G4double E1=minEProj;
-  std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
-  std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
+  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
+  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
   log_ESec_vector->clear();
   log_Prob_vector->clear();
@@ -537 +333 @@
   log_Prob_vector->push_back(-50.);
   
-  G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
+  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
   G4double fE=std::pow(10.,1./nbin_pro_decade);
   G4double int_cross_section=0.;
@@ -544 +340 @@
   
   while (E1 <maxEProj*0.9999999){
-        //G4cout<<E1<<'\t'<<E2<<std::endl;
-        
-        int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
-        //G4cout<<"int_cross_section 1 "<<E1<<'\t'<<int_cross_section<<std::endl;
+        
+        int_cross_section +=integral.Simpson(this,
+        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
         log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
         log_Prob_vector->push_back(std::log(int_cross_section));
@@ -557 +352 @@
   res_mat.clear();
   
-  //if (int_cross_section >0.) {
+  if (int_cross_section >0.) {
         res_mat.push_back(log_ESec_vector);
         res_mat.push_back(log_Prob_vector);
-  //}   
+  }
+  
+ 
   
   return res_mat;
@@ -571 +368 @@
                                 G4double kinEnergyScatProj,
                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
-{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
+{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
   SelectedMaterial= aMaterial;
   kinEnergyScatProjForIntegration = kinEnergyScatProj;
-  /*G4cout<<name<<std::endl;
-  G4cout<<aMaterial->GetName()<<std::endl;
-  G4cout<<kinEnergyScatProj/MeV<<std::endl;*/
+ 
   //compute the vector of integrated cross sections
   //-------------------
@@ -590 +385 @@
   
   
-  std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
-  std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
+  std::vector< double>*  log_ESec_vector = new std::vector< double>();
+  std::vector< double>*  log_Prob_vector = new std::vector< double>();
   log_ESec_vector->push_back(std::log(dEmin));
   log_Prob_vector->push_back(-50.);
-  G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
+  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
   G4double fE=std::pow(dEmax/dEmin,1./nbins);
   
@@ -602 +397 @@
         dE2=dE1*fE;
         int_cross_section +=integral.Simpson(this,
-        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
-        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
-        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
+        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
+        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
         log_Prob_vector->push_back(std::log(int_cross_section));
         dE1=dE2;
@@ -631 +425 @@
   G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
   if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
-  std::vector< G4double >* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
-  //G4double dLog = theMatrix->GetDlog();
-  
-  
+  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
   
   if (theLogPrimEnergyVector->size() ==0){
-        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
-        G4cout<<"The sampling procedure will be stopped."<<std::endl;
+        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
+        G4cout<<"The sampling procedure will be stopped."<<G4endl;
         return 0.;
         
@@ -650 +441 @@
   G4double aLogPrimEnergy1,aLogPrimEnergy2;
   G4double aLogCS1,aLogCS2;
-   G4double log01,log02;
-  std::vector< G4double>* aLogSecondEnergyVector1 =0;
-  std::vector< G4double>* aLogSecondEnergyVector2  =0;
-  std::vector< G4double>* aLogProbVector1=0;
-  std::vector< G4double>* aLogProbVector2=0; 
+  G4double log01,log02;
+  std::vector< double>* aLogSecondEnergyVector1 =0;
+  std::vector< double>* aLogSecondEnergyVector2  =0;
+  std::vector< double>* aLogProbVector1=0;
+  std::vector< double>* aLogProbVector2=0; 
   std::vector< size_t>* aLogProbVectorIndex1=0;
   std::vector< size_t>* aLogProbVectorIndex2=0;
@@ -674 +465 @@
   G4double Emax=0.;
   if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
-        //G4cout<<"Here "<<std::endl;
-        if (ApplyCutInRange) {
+        Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
+        Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
+        G4double dE=0;
+        if (Emin < Emax ){
+            if (ApplyCutInRange) {
                 if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
-                /*if (IsIonisation){
-                        G4double inv_Tcut= 1./currentTcutForDirectSecond;
-                        G4double inv_dE=inv_Tcut-rand_var*(inv_Tcut-1./aPrimEnergy);
-                        Esec= aPrimEnergy+1./inv_dE;
-                        //return Esec;
-                        G4double dE1=currentTcutForDirectSecond;
-                        G4double dE2=currentTcutForDirectSecond*1.00001;
-                        G4double dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
-                        G4double dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
-                        G4double alpha1=std::log(dCS1/dCS2)/std::log(dE1/dE2);
-                        G4double a1=dCS1/std::pow(dE1,alpha1);
-                        dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
-                        dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
-                        
-                        return Esec;
-                        
-                        
-                        
-                        dE1=aPrimEnergy/1.00001;
-                        dE2=aPrimEnergy;
-                        dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
-                        dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
-                        G4double alpha2=std::log(dCS1/dCS2)/std::log(dE1/dE2);
-                        G4double a2=dCS1/std::pow(dE1,alpha1);
-                        return Esec;
-                        
-                        
-                        
-                        
-                }*/     
+                
                 log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
                 log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
                 
-        }       
-        log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
-        log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
-        
-        /*log_dE1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
-        log_dE2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
-        */                                  
-                                            
-        
-        
-        
-        
-        Esec = aPrimEnergy +
-        std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
-        
-        Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy);
-        Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
+            }   
+            log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
+            log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
+             dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
+        }
+        
+        Esec = aPrimEnergy +dE;
         Esec=std::max(Esec,Emin);
         Esec=std::min(Esec,Emax);
         
-        
-        //G4cout<<"Esec "<<Esec<<std::endl;
-        //if (Esec > 2.*aPrimEnergy && second_part_of_same_type) Esec = 2.*aPrimEnergy; 
   }
   else { //Tcut condition is already full-filled
-        /*G4cout<<"Start "<<std::endl;
-        G4cout<<std::exp((*aLogProbVector1)[0])<<std::endl;
-        G4cout<<std::exp((*aLogProbVector2)[0])<<std::endl;*/
-        /*G4double inv_E1= .5/aPrimEnergy;
-                
-        G4double inv_E=inv_E1-rand_var*(inv_E1-0.00001);
-        Esec= 1./inv_E;
-        return Esec;*/
+        
         log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
         log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
-        /*log_E1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
-        log_E2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
-        */
-        
-        
-        /*G4cout<<std::exp(log_E1)<<std::endl;
-        G4cout<<std::exp(log_E2)<<std::endl;*/
         
         Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
@@ -767 +505 @@
  
                                                                                    
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//                              
+G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
+{ SelectCSMatrix(IsScatProjToProjCase);
+  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
+}
+//////////////////////////////////////////////////////////////////////////////
+//
+void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
+{ 
+  indexOfUsedCrossSectionMatrix=0;
+  if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
+  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
+        std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
+        lastCS=lastAdjointCSForScatProjToProjCase;
+        if ( !IsScatProjToProjCase) {
+                CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
+                lastCS=lastAdjointCSForProdToProjCase;
+        }       
+        G4double rand_var= G4UniformRand();
+        G4double SumCS=0.;
+        size_t ind=0;
+        for (size_t i=0;i<CS_Vs_Element->size();i++){
+                SumCS+=(*CS_Vs_Element)[i];
+                if (rand_var<=SumCS/lastCS){
+                        ind=i;
+                        break;
+                }
+        }
+        indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
+  }
 }
 //////////////////////////////////////////////////////////////////////////////
@@ -801 +572 @@
         do {
                 q = G4UniformRand();
-                x = std::pow(xmin, q);
+                x = pow(xmin, q);
                 E=x*Emax;
                 greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1); 
@@ -814 +585 @@
   
   return E;
-  
-  
-  
- 
-                                                                                   
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, 
+                                                  G4double old_weight,  
+                                                  G4double adjointPrimKinEnergy, 
+                                                  G4double projectileKinEnergy,
+                                                  G4bool IsScatProjToProjCase) 
+{
+ G4double new_weight=old_weight;
+ G4double w_corr =1./CS_biasing_factor;
+ w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
+ 
+ 
+ lastCS=lastAdjointCSForScatProjToProjCase;
+ if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
+ if (adjointPrimKinEnergy !=preStepEnergy){ //Is that in all cases needed???
+        G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
+                                                 ,IsScatProjToProjCase );
+        w_corr*=post_stepCS/lastCS;                     
+ }
+        
+ new_weight*=w_corr;
+
+ //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
+ new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
+                                                        //by the factor adjointPrimKinEnergy/projectileKinEnergy
+   
+
+
+ fParticleChange->SetParentWeightByProcess(false);
+ fParticleChange->SetSecondaryWeightByProcess(false);
+ fParticleChange->ProposeParentWeight(new_weight);      
 }
 //////////////////////////////////////////////////////////////////////////////
@@ -830 +630 @@
 //
 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
-{ return PrimAdjEnergy+Tcut;
+{ G4double Emin=PrimAdjEnergy;
+  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
+  return Emin;
 }
 //////////////////////////////////////////////////////////////////////////////
@@ -853 +655 @@
         currentMaterialIndex = currentMaterial->GetIndex();
         size_t idx=56;
-    
+        currentTcutForDirectPrim =0.00000000001;
         if (theAdjEquivOfDirectPrimPartDef) {
                 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
                 else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
                 else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
-                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
-                currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
+                if (idx <56){
+                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
+                        currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
+                }       
         }       
     
-    
+        currentTcutForDirectSecond =0.00000000001;
         if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
                 currentTcutForDirectSecond = currentTcutForDirectPrim;
@@ -873 +677 @@
                         const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
                         currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
+                        if (idx <56){
+                                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
+                                currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
+                        }
+                        
+                        
                 }
         }
   }
 }
+////////////////////////////////////////////////////////////////////////////////////////////
+//
+void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
+{ HighEnergyLimit=aVal;
+  if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
+}
+////////////////////////////////////////////////////////////////////////////////////////////
+//
+void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
+{
+  LowEnergyLimit=aVal;
+  if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
+}
+////////////////////////////////////////////////////////////////////////////////////////////
+//
+void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
+{
+  theAdjEquivOfDirectPrimPartDef=aPart;
+  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
+                                        theDirectPrimaryPartDef=G4Electron::Electron();
+  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
+                                        theDirectPrimaryPartDef=G4Gamma::Gamma();
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4eInverseBremsstrahlung.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 #include "G4eInverseBremsstrahlung.hh"
 #include "G4VEmAdjointModel.hh"
@@ -31 +34 @@
 //
 G4eInverseBremsstrahlung::G4eInverseBremsstrahlung(G4bool whichScatCase,G4String process_name,G4AdjointBremsstrahlungModel* aBremAdjointModel):
-                                G4VAdjointInverseScattering(process_name,whichScatCase)
+                                G4VAdjointReverseReaction(process_name,whichScatCase)
 {theAdjointEMModel = aBremAdjointModel;
- theAdjointEMModel->SetSecondPartOfSameType(false); 
+ theAdjointEMModel->SetSecondPartOfSameType(false);
+ if (IsScatProjToProjCase) SetIntegralMode(true);
+ else   SetIntegralMode(false);
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4eInverseCompton.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 ///////////////////////////////////////////////////////
 // File name:     G4eInverseCompton
@@ -39 +42 @@
 //
 G4eInverseCompton::G4eInverseCompton(G4bool whichScatCase,G4String process_name,G4AdjointComptonModel* aComptonAdjointModel):
-                                G4VAdjointInverseScattering(process_name,whichScatCase)
+                                G4VAdjointReverseReaction(process_name,whichScatCase)
 {theAdjointEMModel = aComptonAdjointModel;
- theAdjointEMModel->SetSecondPartOfSameType(false); 
+ theAdjointEMModel->SetSecondPartOfSameType(false);
+ SetIntegralMode(false);
+ /*if (IsScatProjToProjCase) SetIntegralMode(false);
+ else   SetIntegralMode(true); */
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4eInverseIonisation.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
 ///////////////////////////////////////////////////////
 // File name:     G4eInverseIonisation
@@ -38 +41 @@
 //
 G4eInverseIonisation::G4eInverseIonisation(G4bool whichScatCase,G4String process_name,G4VEmAdjointModel* aEmAdjointModel):
-                                G4VAdjointInverseScattering(process_name,whichScatCase)
+                                G4VAdjointReverseReaction(process_name,whichScatCase)
 {theAdjointEMModel = aEmAdjointModel;
  theAdjointEMModel->SetSecondPartOfSameType(true); 
+ SetIntegralMode(true);
+
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////