00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <fstream>
00038 #include <iomanip>
00039
00040 #include "G4ecpssrFormFactorLixsModel.hh"
00041
00042 #include "globals.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4ios.hh"
00045
00046 #include "G4EMDataSet.hh"
00047 #include "G4LinInterpolation.hh"
00048 #include "G4Proton.hh"
00049 #include "G4Alpha.hh"
00050
00051
00052
00053 G4ecpssrFormFactorLixsModel::G4ecpssrFormFactorLixsModel()
00054 {
00055 interpolation = new G4LinInterpolation();
00056
00057 for (G4int i=6; i<93; i++)
00058 {
00059 protonL1DataSetMap[i] = new G4EMDataSet(i,interpolation);
00060 protonL1DataSetMap[i]->LoadData("pixe/ecpssr/proton/l1-");
00061
00062 protonL2DataSetMap[i] = new G4EMDataSet(i,interpolation);
00063 protonL2DataSetMap[i]->LoadData("pixe/ecpssr/proton/l2-");
00064
00065 protonL3DataSetMap[i] = new G4EMDataSet(i,interpolation);
00066 protonL3DataSetMap[i]->LoadData("pixe/ecpssr/proton/l3-");
00067 }
00068
00069 for (G4int i=6; i<93; i++)
00070 {
00071 alphaL1DataSetMap[i] = new G4EMDataSet(i,interpolation);
00072 alphaL1DataSetMap[i]->LoadData("pixe/ecpssr/alpha/l1-");
00073
00074 alphaL2DataSetMap[i] = new G4EMDataSet(i,interpolation);
00075 alphaL2DataSetMap[i]->LoadData("pixe/ecpssr/alpha/l2-");
00076
00077 alphaL3DataSetMap[i] = new G4EMDataSet(i,interpolation);
00078 alphaL3DataSetMap[i]->LoadData("pixe/ecpssr/alpha/l3-");
00079 }
00080
00081 }
00082
00083
00084
00085 G4ecpssrFormFactorLixsModel::~G4ecpssrFormFactorLixsModel()
00086 {
00087 protonL1DataSetMap.clear();
00088 alphaL1DataSetMap.clear();
00089
00090 protonL2DataSetMap.clear();
00091 alphaL2DataSetMap.clear();
00092
00093 protonL3DataSetMap.clear();
00094 alphaL3DataSetMap.clear();
00095
00096 delete interpolation;
00097 }
00098
00099
00100
00101 G4double G4ecpssrFormFactorLixsModel::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00102 {
00103 G4Proton* aProton = G4Proton::Proton();
00104 G4Alpha* aAlpha = G4Alpha::Alpha();
00105 G4double sigma = 0;
00106
00107 if (energyIncident > 0.1*MeV && energyIncident < 100.*MeV && zTarget < 93 && zTarget > 5) {
00108
00109 if (massIncident == aProton->GetPDGMass())
00110 {
00111 sigma = protonL1DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00112 if (sigma !=0 && energyIncident > protonL1DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00113 }
00114 else if (massIncident == aAlpha->GetPDGMass())
00115 {
00116 sigma = alphaL1DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00117 if (sigma !=0 && energyIncident > alphaL1DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00118 }
00119 else
00120 {
00121 sigma = 0.;
00122 }
00123 }
00124
00125
00126
00127 return sigma;
00128 }
00129
00130
00131
00132 G4double G4ecpssrFormFactorLixsModel::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00133 {
00134 G4Proton* aProton = G4Proton::Proton();
00135 G4Alpha* aAlpha = G4Alpha::Alpha();
00136 G4double sigma = 0;
00137
00138 if (energyIncident > 0.1*MeV && energyIncident < 100.*MeV && zTarget < 93 && zTarget > 5) {
00139
00140 if (massIncident == aProton->GetPDGMass())
00141 {
00142 sigma = protonL2DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00143 if (sigma !=0 && energyIncident > protonL2DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00144 }
00145 else if (massIncident == aAlpha->GetPDGMass())
00146 {
00147 sigma = alphaL2DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00148 if (sigma !=0 && energyIncident > alphaL2DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00149 }
00150 else
00151 {
00152 sigma = 0.;
00153 }
00154 }
00155
00156
00157
00158 return sigma;
00159 }
00160
00161
00162
00163 G4double G4ecpssrFormFactorLixsModel::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00164 {
00165 G4Proton* aProton = G4Proton::Proton();
00166 G4Alpha* aAlpha = G4Alpha::Alpha();
00167 G4double sigma = 0;
00168
00169 if (energyIncident > 0.1*MeV && energyIncident < 100.*MeV && zTarget < 93 && zTarget > 5) {
00170
00171 if (massIncident == aProton->GetPDGMass())
00172 {
00173 sigma = protonL3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00174 if (sigma !=0 && energyIncident > protonL3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00175 }
00176 else if (massIncident == aAlpha->GetPDGMass())
00177 {
00178 sigma = alphaL3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
00179 if (sigma !=0 && energyIncident > alphaL3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
00180 }
00181 else
00182 {
00183 sigma = 0.;
00184 }
00185 }
00186
00187
00188
00189 return sigma;
00190 }