* Change energyloss in NPAnalysis/Gaspard

NPAnalysis/Gaspard/include/ObjectManager.hh

@@ -104,8 +104,9 @@ using namespace NPL ;
 namespace ENERGYLOSS
 {
    // Declare your Energy loss here
-//   EnergyLoss DeutonTargetCD2 = EnergyLoss("deuton_cd2.txt", 100, 1, 2);
-   EnergyLoss DeutonTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 2);
+//   EnergyLoss LightTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 1); // LISE++
+   EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", 100);	// G4
+//   EnergyLoss BeamTarget = EnergyLoss("proton_CD2.G4table", 100);	// G4
 }
 
 using namespace ENERGYLOSS ;

NPAnalysis/Gaspard/src/Analysis.cc

@@ -26,6 +26,9 @@ int main(int argc,char** argv)
    NPL::Reaction* myReaction = new Reaction();
    myReaction->ReadConfigurationFile(reactionfileName);
 
+   // set energy beam at target middle
+   myReaction->SetBeamEnergy(1292);
+
    // Initialize the detector
    NPA::DetectorManager* myDetector = new DetectorManager;
    myDetector->ReadConfigurationFile(detectorfileName);
@@ -81,7 +84,7 @@ int main(int argc,char** argv)
          ThetaStrip = ThetaCalculation (A ,TVector3(0,0,1));
 
          // Correct for energy loss in the target
-         E = DeutonTargetCD2.EvaluateInitialEnergy(E, 5.15*micrometer, ThetaStrip);
+         E = LightTarget.EvaluateInitialEnergy(E, 5.15*micrometer, ThetaStrip);
 
          // Calculate excitation energy
          if (Theta/deg > 90) {

NPSimulation/src/EventGeneratorTransfert.cc

@@ -62,18 +62,20 @@ EventGeneratorTransfert::EventGeneratorTransfert()
    m_ShootLight		= 0;
    m_ShootHeavy		= 0;
 }
+
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void	EventGeneratorTransfert::SetTarget(Target* Target) 
-   {
-   	if(Target!=0)	
-   		{
-   			m_Target = Target;
-   			G4int LightZ = m_Reaction->GetNucleus3()->GetZ() ;
-  			G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
-   			m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZ,LightA, 0.) ,0, m_BeamEnergy);
-   		}
-   
+void EventGeneratorTransfert::SetTarget(Target* Target) 
+{
+   if (Target != 0) {
+      m_Target = Target;
+      G4int LightZ = m_Reaction->GetNucleus3()->GetZ();
+      G4int LightA = m_Reaction->GetNucleus3()->GetA();
+      m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZ,LightA, 0.) ,0, m_BeamEnergy);
    }
+}
+
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 EventGeneratorTransfert::~EventGeneratorTransfert()
 {
@@ -363,16 +365,16 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
    G4double FinalBeamEnergy = 0 ;
    G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
    
-	m_Target->CalculateBeamInteraction(	0, m_SigmaX, 0, m_SigmaThetaX,
-                            					0, m_SigmaY, 0, m_SigmaPhiY,
-				                            	InitialBeamEnergy,
-				                            	BeamName,
-				                           	 	InterCoord, Beam_thetaX, Beam_phiY,
-                            					Beam_theta, Beam_phi,
-				                           	 	FinalBeamEnergy);
+   m_Target->CalculateBeamInteraction(0, m_SigmaX, 0, m_SigmaThetaX,
+                                      0, m_SigmaY, 0, m_SigmaPhiY,
+				      InitialBeamEnergy,
+				      BeamName,
+				      InterCoord, Beam_thetaX, Beam_phiY,
+                                      Beam_theta, Beam_phi,
+				      FinalBeamEnergy);
 				                           	 	
-	m_Reaction->SetBeamEnergy(FinalBeamEnergy);
-  m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
+   m_Reaction->SetBeamEnergy(FinalBeamEnergy);
+   m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
   
    // write vertex position to ROOT file
    G4double x0 = InterCoord.x();

NPSimulation/src/Target.cc

@@ -64,7 +64,7 @@ Target::Target()
    m_WindowsThickness   = 0   ;
    m_TargetTemperature  = 0   ;
    m_TargetPressure  	= 0   ;
-   m_TargetNbLayers      = 50;	// Number of steps by default
+   m_TargetNbLayers     = 50;	// Number of steps by default
 }
 
 G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Temperature, G4double Pressure)
@@ -544,6 +544,8 @@ void Target::CalculateBeamInteraction(	double MeanPosX, double SigmaPosX, double
       // Calculation of effective target thickness and z-position of interaction
       // when the target is tilted wrt the beam axis
       double EffectiveThickness = m_TargetThickness / (BeamDir.unit()).dot(TargetNormal.unit());
+      // remove compilation warning
+      EffectiveThickness *= 1;
       double uniform = RandFlat::shoot();
       z0 = dz + (-m_TargetThickness / 2 + uniform * m_TargetThickness);
 
@@ -557,27 +559,20 @@ void Target::CalculateBeamInteraction(	double MeanPosX, double SigmaPosX, double
       z0 += m_TargetZ;
       InterCoord = G4ThreeVector(x0, y0, z0);
       
-		if(m_TargetType)
-			{
-				G4EmCalculator emCalculator;		
-				if(m_TargetThickness!=0)
-					{
-						for (G4int i = 0; i < m_TargetNbLayers; i++) 
-							{
-								G4double dedx = emCalculator.ComputeTotalDEDX(IncidentBeamEnergy, BeamName, m_TargetMaterial);
-								G4double de   = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers;
-								IncidentBeamEnergy -= de;
-							}
-					
-					}
-			
-			}
-		
-			
-		else
-			{		G4EmCalculator emCalculator;		
-					//	Windows
-					if(m_WindowsThickness!=0)
+      if (m_TargetType) {
+         G4EmCalculator emCalculator;		
+         if (m_TargetThickness != 0) {
+            for (G4int i = 0; i < m_TargetNbLayers; i++) {
+               G4double dedx = emCalculator.ComputeTotalDEDX(IncidentBeamEnergy, BeamName, m_TargetMaterial);
+               G4double de   = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers;
+               IncidentBeamEnergy -= de;
+            }
+         }
+      }			
+      else {
+         G4EmCalculator emCalculator;		
+         // Windows
+         if (m_WindowsThickness != 0)
 					for (G4int i = 0; i < m_TargetNbLayers; i++) 
 						{
 							G4double dedx = emCalculator.ComputeTotalDEDX(IncidentBeamEnergy, BeamName, m_WindowsMaterial);
@@ -596,7 +591,7 @@ void Target::CalculateBeamInteraction(	double MeanPosX, double SigmaPosX, double
 			
 			}
 		
-FinalBeamEnergy=IncidentBeamEnergy;
+   FinalBeamEnergy=IncidentBeamEnergy;
 }
 
 
@@ -621,44 +616,38 @@ void Target::RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double
 }
 
 //	Generate a DEDX file table using the material used in the target
-void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4double Emax)
-	{
-		//	Opening hte output file
-		G4String GlobalPath = getenv("NPTOOL");
-   	G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_TargetMaterial->GetName() + ".G4table";
-	
-		ofstream File		;
-		File.open(Path)	;
-		
-		if(!File) return ;
-		
-		File	<< "Table from Geant4 generate using NPSimulation \t"
-					<< "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_TargetMaterial->GetName() << endl ;
-		
-		G4EmCalculator emCalculator;
+void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin, G4double Emax)
+{
+   // Opening hte output file
+   G4String GlobalPath = getenv("NPTOOL");
+   G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_TargetMaterial->GetName() + ".G4table";
 	
-		for (G4double E=Emin*MeV; E < Emax*MeV; E+=(Emax-Emin)*MeV/10000.) 
-						{
-							G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_TargetMaterial);
-							File << E/MeV << "\t" << dedx/(MeV/um) << endl ;
-						}
-		File.close();
+   ofstream File;
+   File.open(Path);
+
+   if (!File) return;		
+   File << "Table from Geant4 generate using NPSimulation \t"
+	<< "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_TargetMaterial->GetName() << endl;
 		
-		if(!m_TargetType)
-			{
-				G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_WindowsMaterial->GetName() + ".G4table";
-				File.open(Path)		;
-				if(!File) return 	;
-				File	<< "Table from Geant4 generate using NPSimulation \t " 
-					<< "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_WindowsMaterial->GetName() << endl ;
-					
-				for (G4double E=Emin*MeV; E < Emax*MeV; E+=(Emax-Emin)*MeV/10000.) 
-						{
-							G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial);
-							File << E/MeV << "\t" << dedx/(MeV/um) << endl ;
-						}
-			}
-			
-						
-		File.close();
-	}
+   G4EmCalculator emCalculator;	
+   for (G4double E=Emin*MeV; E < Emax*MeV; E+=(Emax-Emin)*MeV/10000.) {
+      G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_TargetMaterial);
+      File << E/MeV << "\t" << dedx/(MeV/um) << endl;
+   }
+   File.close();
+
+   if (!m_TargetType) {
+      G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_WindowsMaterial->GetName() + ".G4table";
+      File.open(Path);
+      if (!File) return;
+      File << "Table from Geant4 generate using NPSimulation \t " 
+           << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_WindowsMaterial->GetName() << endl;
+
+      for (G4double E=Emin*MeV; E < Emax*MeV; E+=(Emax-Emin)*MeV/10000.) {
+         G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial);
+         File << E/MeV << "\t" << dedx/(MeV/um) << endl ;
+      }
+   }
+
+   File.close();
+}