From c4addd3933edd3ed6ff419b3b300e36e5f789cbb Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Thu, 19 Nov 2009 09:39:47 +0000
Subject: [PATCH] * Saving commit, work on energy loss table to be finnish
---
NPAnalysis/10He_Riken/Analysis | Bin 55635 -> 55635 bytes
.../10He_Riken/include/ObjectManager.hh | 56 ++++----
NPAnalysis/10He_Riken/src/Analysis.cc | 2 +-
NPLib/Tools/NPEnergyLoss.cxx | 123 +++++++++---------
.../src/EventGeneratorTransfertToResonance.cc | 5 +
NPSimulation/src/Target.cc | 20 ++-
6 files changed, 104 insertions(+), 102 deletions(-)
diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis
index 75f4cb005797b9560b71bcab202eac2dbc49b5fb..592ab93c584150d90ec51a9f190d63e0e5f2bfa5 100755
GIT binary patch
delta 165
zcmcb-iTUy-<_#5$BDZ%>Zs)9H{&cfFa-Q9ilcD##mn3ZNVO$l)$;iOKaNI?OhhegP
zL>;5e<b4tPNMhfBV&0QWBeyYzP1cFp#+W$yF_0{pJS%z|W7p=un47{19;wFhzLhzN
yMSAWgC5cHnsSF?y7b6&VvQI(a=Gg_hN`f;!uyBB^V`SiG@BlJqZoWCG$_)Tx%sP7j
delta 155
zcmcb-iTUy-<_#5$A_^~y*IkQw!~L$>eO7Xm__s^iAAWD{VO$k9c}b+$WS@vaMvKY&
zBJ{yL=Ey=u-^rzs+ZZDz>qKp1OrHE0NR~{V6}^qIXLDf8O<^TtkJNbI%ACX^y^@L&
z1`yB1C_S+%-Y_ZN$jlogIN7TpaP!;(Z6(23A6PgT85kHC8Tc7IfQ;FjZ%(Rm0{}Y}
BJLdoZ
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index 5657cc8c3..1ab104dbc 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh
@@ -107,12 +107,35 @@ namespace ENERGYLOSS
{
-// // 3He Energy Loss
-// EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" ,
-// 10000 );
+ // 3He Energy Loss
+ EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" ,
+ 10000 );
+
+ EnergyLoss He3TargetGaz = EnergyLoss ( "He3_D2.G4table" ,
+ 10000 );
+
+ EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
+ 10000 ,
+ 1 ,
+ 3 );
+
+ EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
+ 10000 ,
+ 1 ,
+ 3 );
+
+
+
+// // 3He Energy Loss
+// EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
+// 10000 ,
+// 1 ,
+// 3 );
//
-// EnergyLoss He3TargetGaz = EnergyLoss ( "He3_D2.G4table" ,
-// 10000 );
+// EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
+// 10000 ,
+// 1 ,
+// 3 );
//
// EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
// 10000 ,
@@ -124,29 +147,6 @@ namespace ENERGYLOSS
// 1 ,
// 3 );
-
-
- // 3He Energy Loss
- EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
- 10000 ,
- 1 ,
- 3 );
-
- EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
- 10000 ,
- 1 ,
- 3 );
-
- EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
- 10000 ,
- 1 ,
- 3 );
-
- EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
- 10000 ,
- 1 ,
- 3 );
-
// proton Energy Loss
EnergyLoss protonTargetWind = EnergyLoss ( "proton_Mylar.txt" ,
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index 6051ef9c3..ed1a372c4 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -138,7 +138,7 @@ int main(int argc,char** argv)
if(M2 -> GetPositionOfInteraction(hit).Z() > 0)
{
if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0)
- {
+ {
ELab[hit] = M2 -> Si_E[hit] ;
ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
diff --git a/NPLib/Tools/NPEnergyLoss.cxx b/NPLib/Tools/NPEnergyLoss.cxx
index cd3be4860..6285a0247 100644
--- a/NPLib/Tools/NPEnergyLoss.cxx
+++ b/NPLib/Tools/NPEnergyLoss.cxx
@@ -71,58 +71,58 @@ EnergyLoss::EnergyLoss(string Path , int NumberOfSlice=100 , int LiseColumn , i
//If LiseColumn is set to 0 File type is expected to be from SRIM or Geant4
if (LiseColumn == 0)
{
- // Opening dE/dX file
+ // Opening dE/dX file
- ifstream TableFile ;
- TableFile.open(Path.c_str()) ;
+ ifstream TableFile ;
+ TableFile.open(Path.c_str()) ;
- if ( !TableFile )
- {
- cout << "Failed to open file " << Path << endl;
- return;
- }
-
- /* else
- {
- // Reading Data
- double energy, nuclear, electronic;
- string unit, dummy;
-
- while ( TableFile >> energy >> unit
- >> electronic >> nuclear
- >> dummy >> dummy >> dummy
- >> dummy >> dummy >> dummy )
- {
- if ( unit == "keV" ) energy = energy*keV ;
- if ( unit == "MeV" ) energy = energy*MeV ;
- if ( unit == "GeV" ) energy = energy*GeV ;
- fEnergy .push_back ( energy ) ;
- fdEdX_Nuclear .push_back ( nuclear ) ;
- fdEdX_Electronic .push_back ( electronic ) ;
- fdEdX_Total .push_back ( nuclear + electronic ) ;
- }
-
- // Close File
- TableFile.close();
- }*/
-
- else
- {
- // Reading Data
- double energy, total;
- string dummy;
- //skipped first line
- getline(TableFile,dummy);
- while ( TableFile >> energy)
- {
- fEnergy .push_back ( energy*MeV ) ;
- TableFile >> total;
- fdEdX_Total .push_back ( total*MeV/um ) ;
- }
-
- // Close File
- TableFile.close();
- }
+ if ( !TableFile )
+ {
+ cout << "Failed to open file " << Path << endl;
+ return;
+ }
+ else
+ {
+ cout << "Reading Energy Loss File: " << Path << endl ;
+ // Reading Data
+ double energy, total;
+ string dummy;
+ //skipped first line
+ getline(TableFile,dummy);
+ while ( TableFile >> energy)
+ {
+ fEnergy.push_back ( energy*MeV ) ;
+ TableFile >> total;
+ fdEdX_Total.push_back ( total*MeV/micrometer ) ;
+ }
+
+ // Close File
+ TableFile.close();
+ }
+ /* else
+ {
+ // Reading Data
+ double energy, nuclear, electronic;
+ string unit, dummy;
+
+ while ( TableFile >> energy >> unit
+ >> electronic >> nuclear
+ >> dummy >> dummy >> dummy
+ >> dummy >> dummy >> dummy )
+ {
+ if ( unit == "keV" ) energy = energy*keV ;
+ if ( unit == "MeV" ) energy = energy*MeV ;
+ if ( unit == "GeV" ) energy = energy*GeV ;
+ fEnergy .push_back ( energy ) ;
+ fdEdX_Nuclear .push_back ( nuclear ) ;
+ fdEdX_Electronic .push_back ( electronic ) ;
+ fdEdX_Total .push_back ( nuclear + electronic ) ;
+ }
+
+ // Close File
+ TableFile.close();
+ }*/
+
}
@@ -141,7 +141,8 @@ EnergyLoss::EnergyLoss(string Path , int NumberOfSlice=100 , int LiseColumn , i
}
else
- { cout << "Reading Energy Loss File: " << Path << endl ;
+ {
+ cout << "Reading Energy Loss File: " << Path << endl ;
// Reading Data
double energy=0, energyloss=0;
string dummy;
@@ -226,10 +227,8 @@ double EnergyLoss::EvaluateTotalLoss(double Energy) const
return -1000;
}
- Interpolator* s = new Interpolator( fEnergy , fdEdX_Total ) ;
- double val = s->Eval(Energy) ;
-
- delete s ;
+ double val = fInter->Eval(Energy) ;
+
return val ;
}
@@ -267,17 +266,17 @@ double EnergyLoss::Slow( double Energy , // Energy of the detected particle
}
delete s ;
- return slow ;
+ return slow ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-double EnergyLoss::EvaluateInitialEnergy( double Energy , // Energy of the detected particle
- double TargetThickness , // Target Thickness at 0 degree
- double Angle ) // Particle Angle
- const
+double EnergyLoss::EvaluateInitialEnergy( double Energy , // Energy of the detected particle
+ double TargetThickness , // Target Thickness at 0 degree
+ double Angle ) // Particle Angle
+ const
{
// Lise file are given in MeV/u
- // For SRIM file fNumberOfMass = 1 whatever is the nucleus, file are given in MeV
+ // For SRIM and geant4 file fNumberOfMass = 1 whatever is the nucleus, file are given in MeV
Energy = Energy / (double) fNumberOfMass ;
if (Angle > halfpi) Angle = pi-Angle ;
@@ -288,7 +287,7 @@ double EnergyLoss::EvaluateInitialEnergy( double Energy , // Energy of the de
for (int i = 0; i < fNumberOfSlice ; i++)
{
double de = fInter->Eval(Energy) * SliceThickness ;
- Energy += de/fNumberOfMass ;
+ Energy += de/fNumberOfMass ;
}
return (Energy*fNumberOfMass) ;
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index 2317f87b4..06e8d7de2 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -376,6 +376,11 @@ while(ReadingStatus){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
{
+
+ G4int LightZx = m_Reaction->GetNucleus3()->GetZ() ;
+ G4int LightAx = m_Reaction->GetNucleus3()->GetA() ;
+ m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, 550);
+
// Clear contents of Precedent event (now stored in ROOTOutput)
m_InitConditions->Clear();
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index b9e382610..122bbe8bf 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -568,11 +568,8 @@ void Target::CalculateBeamInteraction( double MeanPosX, double SigmaPosX, double
G4double de = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers;
IncidentBeamEnergy -= de;
}
-
}
-
}
-
else
{ G4EmCalculator emCalculator;
@@ -593,7 +590,6 @@ void Target::CalculateBeamInteraction( double MeanPosX, double SigmaPosX, double
G4double de = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers;
IncidentBeamEnergy -= de;
}
-
}
FinalBeamEnergy=IncidentBeamEnergy;
@@ -637,10 +633,10 @@ void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4doub
G4EmCalculator emCalculator;
- for (G4double E=Emin*MeV; E < Emax*MeV; E+=(Emax-Emin)*MeV/10000.)
+ for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/10000.)
{
G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_TargetMaterial);
- File << E/MeV << "\t" << dedx/(MeV/um) << endl ;
+ File << E/MeV << "\t" << dedx/(MeV/micrometer) << endl ;
}
File.close();
@@ -652,13 +648,15 @@ void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4doub
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.)
+ for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/10000.)
{
- G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial);
- File << E/MeV << "\t" << dedx/(MeV/um) << endl ;
+// G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial);
+ G4double dedx = emCalculator.ComputeDEDX( E, Particle ,
+ "ionIoni", m_WindowsMaterial);
+cout << dedx<<endl ;
+ File << E/MeV << "\t" << dedx/(MeV/micrometer) << endl ;
}
}
-
-
File.close();
+
}
--
GitLab