Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointAlongStepWeightCorrection.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointAlongStepWeightCorrection.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointAlongStepWeightCorrection.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,4 +33,7 @@
 #include "G4AdjointCSManager.hh"
 
+#ifdef WIN32
+#include <G4float .h>
+#endif
 
 ///////////////////////////////////////////////////////
@@ -45,16 +51,11 @@
 {; 
 }
-
-
 ///////////////////////////////////////////////////////
 //
-
 void G4AdjointAlongStepWeightCorrection::PreparePhysicsTable(
      const G4ParticleDefinition& )
 {
 ; 
-
 }
-
 ///////////////////////////////////////////////////////
 //
@@ -63,8 +64,4 @@
 {;
 }
-
-
-
-
 ///////////////////////////////////////////////////////
 //
@@ -85,7 +82,26 @@
   									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,2 +112,12 @@
 
 }
+///////////////////////////////////////////////////////
+//
+G4double G4AdjointAlongStepWeightCorrection::GetContinuousStepLimit(const G4Track& track,
+                G4double , G4double , G4double& )
+{ 
+  G4double x = DBL_MAX;
+  DefineMaterial(track.GetMaterialCutsCouple());
+  preStepKinEnergy = track.GetKineticEnergy();
+  return x;
+}
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,197 +33,265 @@
 #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,194 +304,57 @@
  
  
+ //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,237 +362,33 @@
 ////////////////////////////////////////////////////////////////////////////////
 //
-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;	
+}					     
+
+
+
+
+
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,8 +43,12 @@
 #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,16 +72,30 @@
   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,5 +117,5 @@
   if (anAdjPartDef && aProcess){
   	RegisterAdjointParticle(anAdjPartDef);
-	int index=-1;
+	G4int index=-1;
 	
   	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
@@ -111,5 +132,5 @@
   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,5 +142,5 @@
 //
 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,4 +153,9 @@
 	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,7 +188,9 @@
 	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,5 +201,5 @@
 				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,8 +208,8 @@
 					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,5 +221,5 @@
 					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,5 +231,6 @@
 			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,6 +238,10 @@
 		}		
 	}
-	G4cout<<"All adjoint cross section matrices are built "<<std::endl;
+	G4cout<<"              All adjoint cross section matrices are computed!"<<G4endl;
+	G4cout<<"======================================================================"<<G4endl;
+	
 	CrossSectionMatrixesAreBuilt = true;
+
+
 }
 
@@ -221,13 +254,40 @@
   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,20 +295,63 @@
 			}
 			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,12 +365,7 @@
 	 		   	     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,13 +377,85 @@
 	 		   	     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,54 +463,54 @@
 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,5 +551,5 @@
 	CS_Vs_Element.clear();
 	if (!aModel->GetUseMatrix()){
-		return aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase);
+		CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
 				         
 	
@@ -397,16 +567,7 @@
 					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,5 +578,5 @@
 			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,5 +587,5 @@
 				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,16 +614,8 @@
   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,5 +626,5 @@
 						        	   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,28 +651,11 @@
 {
  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,10 +663,12 @@
 		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,9 +677,7 @@
  	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,6 +685,4 @@
 					    	       listOfAdjointEMModel[i], 
 					    	       Ekin, Tlow,false, CS_Vs_Element);
-						       
-			//G4cout<<"Yes2 "<<TotalCS<<std::endl;
 		}
 		
@@ -563,5 +697,5 @@
 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,34 +711,23 @@
    
     
-   
-   
-   
-   
-   
    //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,22 +739,20 @@
    //-----------------------------------------
    
-   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,17 +762,11 @@
    
    
-   
-   
-   
-   
-   
   std::vector<G4AdjointCSMatrix*> res;
   res.clear();
-  
   res.push_back(theCSMatForProdToProjBackwardScattering);
   res.push_back(theCSMatForScatProjToProjBackwardScattering);
   
 
-#ifdef TEST_MODE
+/*
   G4String file_name;
   std::stringstream astream;
@@ -662,12 +777,6 @@
   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,23 +811,23 @@
    //-----------------------------------------
    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,22 +842,23 @@
    //-----------------------------------------
    
-   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,8 +880,8 @@
   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,5 +896,8 @@
 {
  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,5 +907,7 @@
 {
  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,23 +920,41 @@
     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,15 +967,11 @@
  
   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,16 +980,8 @@
   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;
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc	(revision 1196)
@@ -24,10 +24,12 @@
 // ********************************************************************
 //
+// $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,20 +66,18 @@
 ///////////////////////////////////////////////////////
 //
- 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,4 +85,5 @@
 		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,20 +103,14 @@
 ///////////////////////////////////////////////////////
 //
-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,9 +122,9 @@
  	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,11 +132,11 @@
 			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,9 +144,9 @@
 			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,6 +172,6 @@
 		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++){
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,4 +37,5 @@
 #include "G4AdjointGamma.hh"
 #include "G4Gamma.hh"
+#include "G4KleinNishinaCompton.hh"
 
 
@@ -44,5 +48,4 @@
   SetUseMatrix(true);
   SetUseMatrixPerElement(true);
-  SetIsIonisation(false);
   SetUseOnlyOneMatrixForAllElements(true);
   theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
@@ -50,4 +53,5 @@
   theDirectPrimaryPartDef=G4Gamma::Gamma();
   second_part_of_same_type=false;
+  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
  
 }
@@ -62,5 +66,5 @@
 	               G4ParticleChange* fParticleChange)
 { 
-   
+   if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
    
    //A recall of the compton scattering law is 
@@ -71,25 +75,14 @@
    
   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,5 +101,5 @@
  //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,5 +109,5 @@
  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,15 +129,19 @@
  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,5 +151,125 @@
  
  	
-} 			
+} 
+////////////////////////////////////////////////////////////////////////////////
+//
+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,7 +296,9 @@
                                       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,8 +310,8 @@
  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,12 +341,12 @@
  //-------------------------------
  
- 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,4 +354,5 @@
 
 }
+				      
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -251,2 +371,38 @@
   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;	
+}
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,5 +64,5 @@
 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,8 +72,8 @@
 { 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,5 +101,5 @@
   }
   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,7 +107,7 @@
 ///////////////////////////////////////////////////////
 //
-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,8 +152,8 @@
 ///////////////////////////////////////////////////////
 //
-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,9 +173,9 @@
 ///////////////////////////////////////////////////////
 //
-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,5 +183,5 @@
 ///////////////////////////////////////////////////////
 //
-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,14 +205,14 @@
 ///////////////////////////////////////////////////////
 //
-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]);
Index: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointPhotoElectricModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4AdjointPhotoElectricModel.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4AdjointPhotoElectricModel.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,4 +38,5 @@
 #include "G4AdjointGamma.hh"
 
+
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -41,6 +45,13 @@
 
 { 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,37 +68,22 @@
 
   //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,6 +92,6 @@
    //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,5 +109,5 @@
    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,7 +127,7 @@
   
      
-  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,5 +137,5 @@
   //Weight correction
  //-----------------------					   
-  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy);	
+  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase);	
  
   
@@ -149,11 +145,34 @@
   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,10 +183,15 @@
 				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,7 +200,14 @@
 		xsec[index_element] = totAdjointCS;
 	} 
-  }	
-  return totAdjointCS;
-  
+
+	totBiasedAdjointCS=std::min(totAdjointCS,0.01);
+//	totBiasedAdjointCS=totAdjointCS;
+	factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
+	lastCS=totBiasedAdjointCS;
+  	
+	
+  }
+  return totBiasedAdjointCS;
+
   
 }
@@ -188,23 +219,26 @@
   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;
   
Index: trunk/source/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,4 +34,6 @@
 #include "G4VParticleChange.hh"
 #include "G4UnitsTable.hh"
+#include "G4AdjointCSManager.hh"
+#include "G4LossTableManager.hh"
 
 
@@ -38,4 +43,5 @@
   G4ProcessType type): G4VContinuousProcess(name, type)
 {
+
 
   linLossLimit=0.05;
@@ -44,4 +50,14 @@
   is_integral = false;
   
+  //Will be properly set in SetDirectParticle()
+  IsIon=false;
+  massRatio =1.;
+  chargeSqRatio=1.;
+  preStepChargeSqRatio=1.;
+
+  
+
+  
+  
 }
 
@@ -70,6 +86,18 @@
 }
 
-
-
+///////////////////////////////////////////////////////
+//
+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,4 +108,11 @@
 {
    
+  //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,15 +121,12 @@
   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,50 +135,96 @@
   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,5 +232,5 @@
   //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,2 +246,66 @@
   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);
+}
Index: trunk/source/processes/electromagnetic/adjoint/src/G4InversePEEffect.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4InversePEEffect.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4InversePEEffect.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,7 +34,10 @@
 //
 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);
+  
+ 
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Index: trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc	(revision 1196)
@@ -24,8 +24,9 @@
 // ********************************************************************
 //
+// $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,4 +34,5 @@
 #include "G4AdjointElectron.hh"
 #include "G4AdjointInterpolator.hh"
+#include "G4PhysicsTable.hh"
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -39,14 +41,8 @@
 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,138 +50,4 @@
 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,34 +57,19 @@
 { 
   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,5 +77,5 @@
 ////////////////////////////////////////////////////////////////////////////////
 //
-//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,9 +92,6 @@
 	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,18 +100,5 @@
  	
 	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,31 +138,11 @@
 	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,51 +167,37 @@
 //
 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,7 +208,8 @@
                                 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,6 +220,6 @@
   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,5 +227,5 @@
   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,8 +234,8 @@
   
   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,7 +261,7 @@
                                 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,10 +277,14 @@
   
   
-  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,15 +293,11 @@
   	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,10 +317,8 @@
 				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,6 +326,6 @@
   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,5 +333,5 @@
   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,8 +340,7 @@
   
   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,8 +352,10 @@
   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,10 +368,8 @@
 				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,9 +385,9 @@
   
   
-  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,7 +397,6 @@
   	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,12 +425,9 @@
   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,9 +441,9 @@
   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,82 +465,29 @@
   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,4 +505,37 @@
  
 										   
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//				
+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,5 +572,5 @@
 	do {
         	q = G4UniformRand();
-        	x = std::pow(xmin, q);
+        	x = pow(xmin, q);
 		E=x*Emax;
 		greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1); 
@@ -814,9 +585,38 @@
   
   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,5 +630,7 @@
 //
 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
-{ return PrimAdjEnergy+Tcut;
+{ G4double Emin=PrimAdjEnergy;
+  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
+  return Emin;
 }
 //////////////////////////////////////////////////////////////////////////////
@@ -853,14 +655,16 @@
     	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,6 +677,35 @@
     			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();
+}
Index: trunk/source/processes/electromagnetic/adjoint/src/G4eInverseBremsstrahlung.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4eInverseBremsstrahlung.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4eInverseBremsstrahlung.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,7 +34,9 @@
 //
 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);
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Index: trunk/source/processes/electromagnetic/adjoint/src/G4eInverseCompton.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4eInverseCompton.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4eInverseCompton.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,7 +42,10 @@
 //
 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); */
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Index: trunk/source/processes/electromagnetic/adjoint/src/G4eInverseIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/adjoint/src/G4eInverseIonisation.cc	(revision 966)
+++ trunk/source/processes/electromagnetic/adjoint/src/G4eInverseIonisation.cc	(revision 1196)
@@ -24,4 +24,7 @@
 // ********************************************************************
 //
+// $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,7 +41,9 @@
 //
 G4eInverseIonisation::G4eInverseIonisation(G4bool whichScatCase,G4String process_name,G4VEmAdjointModel* aEmAdjointModel):
-				G4VAdjointInverseScattering(process_name,whichScatCase)
+				G4VAdjointReverseReaction(process_name,whichScatCase)
 {theAdjointEMModel = aEmAdjointModel;
  theAdjointEMModel->SetSecondPartOfSameType(true); 
+ SetIntegralMode(true);
+
 }
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
