diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh index 5657cc8c3fdea74b3af5aaa02515548b9eeb5ec0..336fb8c7142848ae3f053a591911a9684948e3f2 100644 --- a/NPAnalysis/10He_Riken/include/ObjectManager.hh +++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh @@ -107,12 +107,39 @@ namespace ENERGYLOSS { -// // 3He Energy Loss -// EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" , -// 10000 ); + // 3He Energy Loss + EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" , + "G4Table", + 10000 ); + + EnergyLoss He3TargetGaz = EnergyLoss ( "He3_D2.G4table" , + "G4Table", + 10000 ); + + EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" , + "LISE", + 10000 , + 1 , + 3 ); + + EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" , + "LISE", + 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,42 +151,22 @@ 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" , + "LISE", 1000 , 1 , 1 ); EnergyLoss protonTargetGaz = EnergyLoss ( "proton_D2gaz_1b_26K.txt" , + "LISE", 1000 , 1 , 1 ); EnergyLoss protonStripAl = EnergyLoss ( "proton_Al.txt" , + "LISE", 100 , 1 , 1 ); diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc index 6051ef9c3be91100f033786a0192e4edd8c37ab3..ed1a372c4b79fa6f9a43e552785d8beb9bd4eb98 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 cd3be4860ae1c5c0f9c69a5e832014dd32bf9c57..f41ac75f9b09e671ac9c8875d142358d2934c5a2 100644 --- a/NPLib/Tools/NPEnergyLoss.cxx +++ b/NPLib/Tools/NPEnergyLoss.cxx @@ -56,7 +56,7 @@ EnergyLoss::~EnergyLoss() {} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -EnergyLoss::EnergyLoss(string Path , int NumberOfSlice=100 , int LiseColumn , int NumberOfMass) +EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int LiseColumn , int NumberOfMass) { fNumberOfSlice = NumberOfSlice ; @@ -67,23 +67,35 @@ EnergyLoss::EnergyLoss(string Path , int NumberOfSlice=100 , int LiseColumn , i cout << "///////////////////////////////// " << endl ; cout << "Initialising an EnergyLoss object " << endl ; - - //If LiseColumn is set to 0 File type is expected to be from SRIM or Geant4 - if (LiseColumn == 0) - { - // Opening dE/dX file - + ifstream TableFile ; TableFile.open(Path.c_str()) ; - if ( !TableFile ) - { - cout << "Failed to open file " << Path << endl; - return; - } - - /* else - { + // Opening dE/dX file + if(!TableFile) cout << "ERROR: TABLE FILE NOT FOUND" << endl; + + if (Source == "G4Table") + { + 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 if (Source == "SRIM") + { // Reading Data double energy, nuclear, electronic; string unit, dummy; @@ -104,67 +116,37 @@ EnergyLoss::EnergyLoss(string Path , int NumberOfSlice=100 , int LiseColumn , i // Close File TableFile.close(); - }*/ - - else - { + } + + else if(Source == "LISE") + { + cout << "Reading Energy Loss File: " << Path << endl ; // Reading Data - double energy, total; - string dummy; - //skipped first line + double energy=0, energyloss=0; + string dummy; + // skipping comment first line getline(TableFile,dummy); - while ( TableFile >> energy) - { - fEnergy .push_back ( energy*MeV ) ; - TableFile >> total; - fdEdX_Total .push_back ( total*MeV/um ) ; + + while ( TableFile >> energy ) + { + for (int k = 0 ; k < 11 ; k++ ) + { + TableFile >> dummy ; + if (k+1==LiseColumn) energyloss = atof(dummy.c_str()) ; + } + fEnergy.push_back (energy*MeV) ; + fdEdX_Total.push_back(energyloss*MeV/micrometer); } // Close File TableFile.close(); - } + } - } - - //Else File is expected to be from Lise, and LiseColumn gives which model to take else { - // Opening dE/dX file - - ifstream TableFile ; - TableFile.open(Path.c_str()) ; - - if ( !TableFile ) - { - cout << "Failed to open file " << Path << endl; - return; - } - - else - { cout << "Reading Energy Loss File: " << Path << endl ; - // Reading Data - double energy=0, energyloss=0; - string dummy; - // skipping comment first line - getline(TableFile,dummy); - - while ( TableFile >> energy ) - { - for (int k = 0 ; k < 11 ; k++ ) - { - TableFile >> dummy ; - if (k+1==LiseColumn) energyloss = atof(dummy.c_str()) ; - } - fEnergy.push_back (energy*MeV) ; - fdEdX_Total.push_back(energyloss*MeV/micrometer); - } - - // Close File - TableFile.close(); - } - - + cout << "ERROR : Wrong Source Type" << endl ; } + fInter = new Interpolator( fEnergy , fdEdX_Total ) ; cout << "///////////////////////////////// " << endl ; } @@ -226,10 +208,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 +247,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 +268,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/NPLib/Tools/NPEnergyLoss.h b/NPLib/Tools/NPEnergyLoss.h index 005b0b2b8ce447f763edf9814a1a25555be16c46..080c3298a0ccd631bff74667bd76eda58f40b028 100644 --- a/NPLib/Tools/NPEnergyLoss.h +++ b/NPLib/Tools/NPEnergyLoss.h @@ -60,10 +60,11 @@ namespace NPL public : // Constructor EnergyLoss(); - EnergyLoss( string Path , // Path of dE/dX table file - int NumberOfSlice , // Low number = Faster, High Number = more accurate / typical: 100 to 1000 - int LiseColumns=0 , // Indicate which model to read in a lise File, set to 0 (Default value) for a SRIM file - int NumberOfMass=1 ); // Number of mass A of the nucleus (used only for Lise file) + EnergyLoss( string Path , // Path of dE/dX table file + string Source , // Type of file : Geant4,Lise,SRIM + int NumberOfSlice , // Low number = Faster, High Number = more accurate / typical: 100 to 1000 + int LiseColumns=0 , // Indicate which model to read in a lise File, set to 0 (Default value) for a SRIM file + int NumberOfMass=1 );// Number of mass A of the nucleus (used only for Lise file) ~EnergyLoss(); private : // dE/dX, slice parameter diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc index b3da7f51cfae0e30bb09efac0fb64eab9350296d..a8c8db2784d1220bd6c98a6e5744e76339d8e86f 100644 --- a/NPSimulation/src/EventGeneratorBeam.cc +++ b/NPSimulation/src/EventGeneratorBeam.cc @@ -185,6 +185,11 @@ void EventGeneratorBeam::ReadConfiguration(string Path) //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun) { + //--------------write the DeDx Table ------------------- + if(m_Target!=0) + m_Target->WriteDEDXTable(m_particle ,0, m_BeamEnergy+4*m_BeamEnergySpread); + + m_InitConditions->Clear(); /////////////////////////////////////////////////////////////////////// diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc index 0a704a9c75f62badacb8846643892616c92a84bd..e30585f7f0bace34e7f3eeea56fda8ba6e05d8da 100644 --- a/NPSimulation/src/EventGeneratorTransfert.cc +++ b/NPSimulation/src/EventGeneratorTransfert.cc @@ -328,6 +328,27 @@ while(ReadingStatus){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun) { + // If first time, write the DeDx table + if(anEvent->GetEventID()==0) + { + //-------------- Before living, wrtie the DeDx Table ------------------- + + G4int LightZx = m_Reaction->GetNucleus3()->GetZ() ; + G4int LightAx = m_Reaction->GetNucleus3()->GetA() ; + + G4int BeamZx = m_Reaction->GetNucleus1()->GetZ() ; + G4int BeamAx = m_Reaction->GetNucleus1()->GetA() ; + + if(m_Target!=0) + { + m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread); + m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(BeamZx,BeamAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread); + } + + } + + + // Clear contents of Precedent event (now stored in ROOTOutput) m_InitConditions->Clear(); diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc index 2317f87b46dd7c557bf08742b6c3ca504cff7701..270f2eadae742b4d751c5211cbee822406692ce2 100644 --- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc +++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc @@ -69,7 +69,7 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance() //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance() { - //------------- Default Destructor ------------ + //------------- Default Destructor ------------ delete m_InitConditions; delete m_Reaction ; } @@ -77,13 +77,7 @@ EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance() void EventGeneratorTransfertToResonance::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...... EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string name1 , @@ -376,6 +370,26 @@ while(ReadingStatus){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun) { + + // If first time, write the DeDx table + if(anEvent->GetEventID()==0) + { + //-------------- Before living, wrtie the DeDx Table ------------------- + + G4int LightZx = m_Reaction->GetNucleus3()->GetZ() ; + G4int LightAx = m_Reaction->GetNucleus3()->GetA() ; + + G4int BeamZx = m_Reaction->GetNucleus1()->GetZ() ; + G4int BeamAx = m_Reaction->GetNucleus1()->GetA() ; + + if(m_Target!=0) + { + m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread); + m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(BeamZx,BeamAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread); + } + + } + // 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 948505403ab45e764c4be906b1553448747aa562..51eba32448deaa208974316e2729fddc5e7b8902 100644 --- a/NPSimulation/src/Target.cc +++ b/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) @@ -292,7 +292,7 @@ void Target::ReadConfiguration(string Path) cout << m_TargetZ / mm << " )" << endl ; } - else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) { + else if (DataBuffer.compare(0, 9, "NbLayers=") == 0) { check_m_TargetNbLayers = true ; ConfigFile >> DataBuffer; m_TargetNbLayers = atoi(DataBuffer.c_str()); @@ -386,7 +386,7 @@ void Target::ReadConfiguration(string Path) cout << m_TargetZ / mm << " )" << endl ; } - else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) { + else if (DataBuffer.compare(0, 9, "m_TargetNbLayers=") == 0) { check_m_TargetNbLayers = true ; ConfigFile >> DataBuffer; m_TargetNbLayers = atoi(DataBuffer.c_str()); @@ -544,8 +544,6 @@ 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); @@ -559,20 +557,24 @@ 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); @@ -588,10 +590,9 @@ void Target::CalculateBeamInteraction( double MeanPosX, double SigmaPosX, double G4double de = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers; IncidentBeamEnergy -= de; } - } - FinalBeamEnergy=IncidentBeamEnergy; +FinalBeamEnergy=IncidentBeamEnergy; } @@ -616,38 +617,45 @@ 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"; +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; + ofstream File ; + File.open(Path) ; - 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(); -} + if(!File) return ; + + File << "Table from Geant4 generate using NPSimulation \t" + << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_TargetMaterial->GetName() << endl ; + + G4EmCalculator emCalculator; + + for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/10000.) + { + G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_TargetMaterial); + File << E/MeV << "\t" << dedx/(MeV/micrometer) << 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; E < Emax; E+=(Emax-Emin)/10000.) + { +// G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial); + G4double dedx = emCalculator.ComputeDEDX( E, Particle , + "ionIoni", m_WindowsMaterial); + File << E/MeV << "\t" << dedx/(MeV/micrometer) << endl ; + } + } + File.close(); + + }