#include "SteppingAction.hh" #include "EventAction.hh" #include "TrackingAction.hh" #include "Transport.hh" #include "G4Step.hh" #include "G4Event.hh" #include "G4RunManager.hh" #include "G4LogicalVolume.hh" #include "G4VPhysicalVolume.hh" #include "G4UnitsTable.hh" #include "G4SystemOfUnits.hh" #include "G4Box.hh" #include "HistoManager.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... SteppingAction::SteppingAction(EventAction* eventAction, TrackingAction* trackingAction) : G4UserSteppingAction(), fEventAction(eventAction), fTrackingAction(trackingAction), fScoringVolume(0) { radioactiveTrackIDEpi1 = -1; radioactiveTrackIDEpi2 = -1; EpiThick = -1; SubThick = -1; histoManager = HistoManager::GetHistoManager(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... SteppingAction::~SteppingAction() {} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void SteppingAction::UserSteppingAction(const G4Step* step) { ///G4cout << "Step in volume " << step->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " " << step->GetPreStepPoint()->GetPosition() << G4endl; // collect energy deposited in this step G4double edepStep = step->GetTotalEnergyDeposit(); fEventAction->AddEdep(edepStep); G4int evtNumber = 0; const G4Event* evt = G4RunManager::GetRunManager()->GetCurrentEvent(); if(evt) evtNumber = evt->GetEventID(); G4double rndmSeed = G4Random::getTheSeed(); G4double edep = step->GetTotalEnergyDeposit(); G4double niel = step->GetNonIonizingEnergyDeposit(); G4StepPoint* prePoint = step->GetPreStepPoint(); G4ThreeVector prePointPosition = prePoint->GetPosition(); G4VPhysicalVolume* pv = prePoint->GetPhysicalVolume(); G4String pvname = pv->GetName(); bool isIncident = false; if( prePoint->GetStepStatus() == fGeomBoundary ){ isIncident = true;} G4StepPoint* postPoint = step->GetPostStepPoint(); G4ThreeVector postPointPosition = postPoint->GetPosition(); bool isExit = false; if( postPoint->GetStepStatus() == fGeomBoundary ) isExit = true; //G4String processName = postPoint->GetProcessDefinedStep()->GetProcessName(); //G4int processType = GetProcessType(processName); G4int processSubType = postPoint->GetProcessDefinedStep()->GetProcessSubType(); G4Track* track = step->GetTrack(); G4int trackID = track->GetTrackID(); G4int parentID = track->GetParentID(); //std::cout << trackID << " " << processSubType << std::endl; const G4DynamicParticle* part = track->GetDynamicParticle(); G4int pdgid = part->GetPDGcode(); G4int charge = part->GetCharge(); //G4double kineticEnergy = part->GetKineticEnergy(); G4double kineticEnergy = prePoint->GetKineticEnergy(); G4ThreeVector momentum = part->GetMomentumDirection(); G4String definition = part->GetDefinition()->GetParticleName(); G4double x1 = prePointPosition.x(); G4double x2 = postPointPosition.x(); G4double x = x1 + G4UniformRand()*(x2-x1); G4double y1 = prePointPosition.y(); G4double y2 = postPointPosition.y(); G4double y = y1 + G4UniformRand()*(y2-y1); G4double z1 = prePointPosition.z(); G4double z2 = postPointPosition.z(); G4double z = z1 + G4UniformRand()*(z2-z1); // kill all tracks that are past 3m to speed up simulation // NOTE: Only valid as dont care about back scatter from walls at the moment if(z/mm > 2000){ track->SetTrackStatus(fStopAndKill); } if(step->GetPreStepPoint()->GetMomentumDirection().z() < 0){ track->SetTrackStatus(fStopAndKill); } if(pvname == "JonMonitor" && histoManager){ histoManager->FillTree("jonMonitor", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } /*if(pvname == "Stiffener_v4" && histoManager){ histoManager->FillTree(pvname, rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); }*/ if(pvname == "ASICs" && histoManager) { G4TouchableHistory*theTouchable = (G4TouchableHistory*)(prePoint->GetTouchable()); G4int motherCopyNo = theTouchable->GetReplicaNumber(1); // get depth==1 to get mother replica number in //G4cout << pvname << " Copy = " << motherCopyNo << " z = " << z << G4endl; G4String strip_names[12] = {"ASICs_x1", "ASICs_u1", "ASICs_v1", "ASICs_x2", "ASICs_u2", "ASICs_v2","ASICs_x3", "ASICs_u3", "ASICs_v3","ASICs_x4", "ASICs_u4", "ASICs_v4"}; if(motherCopyNo>5) motherCopyNo -= 4; //G4cout << motherCopyNo << G4endl; //G4cout << strip_names[motherCopyNo] << G4endl; if(isIncident){ histoManager->FillHisto(strip_names[motherCopyNo]+"_kineticEnergy", kineticEnergy); if(pdgid==2212){ histoManager->FillHisto(strip_names[motherCopyNo]+"_kineticEnergy_proton", kineticEnergy); } G4double neq(0); if(pdgid==2212) neq = GetProton1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==2112) neq = GetNeutron1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==-11 || pdgid==11) neq = GetElectron1MeVNeutronEquivalence(kineticEnergy/MeV); if(neq>0){ histoManager->FillHisto2D(strip_names[motherCopyNo]+"_1MeVNeutronEquivalent", x, y, neq); if(pdgid==2212){ histoManager->FillHisto2D(strip_names[motherCopyNo]+"_1MeVNeutronEquivalent_proton", x, y, neq); } } } histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Edep", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Inelastic", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Elastic", x/mm,y/mm,1); if(pdgid==2212){ histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Edep_proton", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Inelastic_proton", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D(strip_names[motherCopyNo]+"_Elastic_proton", x/mm,y/mm,1); } /* if(motherCopyNo==0){ histoManager->FillTree("ASICs_x1", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==1){ histoManager->FillTree("ASICs_u1", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==2){ histoManager->FillTree("ASICs_v1", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==3){ histoManager->FillTree("ASICs_x2", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==4){ histoManager->FillTree("ASICs_u2", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==5){ histoManager->FillTree("ASICs_v2", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==10){ histoManager->FillTree("ASICs_x3", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==11){ histoManager->FillTree("ASICs_u3", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==12){ histoManager->FillTree("ASICs_v3", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==13){ histoManager->FillTree("ASICs_x4", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==14){ histoManager->FillTree("ASICs_u4", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(motherCopyNo==15){ histoManager->FillTree("ASICs_v4", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); }*/ } if(pvname == "DummySi" && histoManager) { G4TouchableHistory*theTouchable = (G4TouchableHistory*)(prePoint->GetTouchable()); G4int motherCopyNo = theTouchable->GetReplicaNumber(1); // get depth==1 to get mother replica number G4String layerNo = G4UIcommand::ConvertToString(motherCopyNo); if(isIncident){ histoManager->FillHisto("DummySi_RT"+layerNo+"_kineticEnergy", kineticEnergy); if(pdgid==2212){ histoManager->FillHisto("DummySi_RT"+layerNo+"_kineticEnergy_proton", kineticEnergy); } G4double neq(0); if(pdgid==2212) neq = GetProton1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==2112) neq = GetNeutron1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==-11 || pdgid==11) neq = GetElectron1MeVNeutronEquivalence(kineticEnergy/MeV); if(neq>0) histoManager->FillHisto2D("DummySi_RT"+layerNo+"_1MeVNeutronEquivalent", x, y, neq); if(pdgid==2212){ histoManager->FillHisto2D("DummySi_RT"+layerNo+"_1MeVNeutronEquivalent_proton", x, y, neq); } } histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Edep", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Inelastic", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Elastic", x/mm,y/mm,1); if(pdgid==2212){ histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Edep_proton", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Inelastic_proton", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D("DummySi_RT"+layerNo+"_Elastic_proton", x/mm,y/mm,1); } /*histoManager->FillTree("DummySi_RT"+layerNo, rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); */ // look at just the electronics at end of long face if( (x/mm < -65.2*mm && x/mm > -70.0*mm) || (x/mm > 65.2*mm && x/mm < 70.0*mm) ){ if(isIncident){ histoManager->FillHisto("Components_RT"+layerNo+"_kineticEnergy", kineticEnergy); if(pdgid==2212){ histoManager->FillHisto("Components_RT"+layerNo+"_kineticEnergy_proton", kineticEnergy); } G4double neq(0); if(pdgid==2212) neq = GetProton1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==2112) neq = GetNeutron1MeVNeutronEquivalence(kineticEnergy/MeV); if(pdgid==-11 || pdgid==11) neq = GetElectron1MeVNeutronEquivalence(kineticEnergy/MeV); if(neq>0) histoManager->FillHisto2D("Components_RT"+layerNo+"_1MeVNeutronEquivalent", x, y, neq); if(pdgid==2212){ histoManager->FillHisto2D("Components_RT"+layerNo+"_1MeVNeutronEquivalent_proton", x, y, neq ); } } histoManager->FillHisto2D("Components_RT"+layerNo+"_Edep", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D("Components_RT"+layerNo+"_Inelastic", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D("Components_RT"+layerNo+"_Elastic", x/mm,y/mm,1); if(pdgid==2212){ histoManager->FillHisto2D("Components_RT"+layerNo+"_Edep_proton", x/mm,y/mm,edep/MeV); if(processSubType==121) histoManager->FillHisto2D("Components_RT"+layerNo+"_Inelastic_proton", x/mm,y/mm,1); if(processSubType==111) histoManager->FillHisto2D("Components_RT"+layerNo+"_Elastic_proton", x/mm,y/mm,1); } } } /* if(pvname == "Phantom" && histoManager){ histoManager->FillTree("Phantom", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(pvname == "PhantomCompensator" && histoManager){ histoManager->FillTree("PhantomCompensator", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } */ if(pvname.contains("TP1") && isIncident && histoManager){ histoManager->FillTree("TruthPlane_I", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } if(pvname.contains("TP4") && isIncident && histoManager){ histoManager->FillTree("TruthPlane_O", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } // /// // /// TEST to look at the radioactive decay locations // /// /* // if not a primary particle G4int creatorProcess = -1; G4String pname = ""; if(track->GetCreatorProcess()!=NULL) { creatorProcess = track->GetCreatorProcess()->GetProcessSubType(); pname = track->GetCreatorProcess()->GetProcessName(); } if (pname == "RadioactiveDecay" && histoManager) { //G4cout <<"RADIOACTIVE DECAY!!! " << processSubType <FillTree("radioactiveDecay", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, creatorProcess); } if (pname == "Decay" && histoManager) { G4cout <<"DECAY!!! "<FillTree("Decay", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, processSubType); } /////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////// // // // NOTE: This section of the code looks at the hits in the epi which arise from // // radioactive decays. It then stores all of the hits within the DD mother // // that have a creator process (so not a primary). The source of the hits // // in the DD system can then be calculated using this information // // // /////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////// if(pvname == "ddEpi1" && histoManager) { histoManager->FillTree("ddEpi1", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, creatorProcess); if(creatorProcess!=-1) { histoManager->FillTree("radioactiveDecay_ddEpi1", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, creatorProcess); G4cout << "evt number = " << evtNumber << G4endl; G4cout << track->GetPosition() << G4endl; G4cout << track->GetVertexPosition() << G4endl << G4endl; if(radioactiveTrackIDEpi1!=trackID) // only want to count once the source of the decay { radioactiveTrackIDEpi1 = trackID; histoManager->FillTree("radioactiveDecay_ddEpi1_Vertex", rndmSeed, evtNumber, track->GetVertexPosition()[0],track->GetVertexPosition()[1],track->GetVertexPosition()[2], track->GetVertexMomentumDirection()[0],track->GetVertexMomentumDirection()[1],track->GetVertexMomentumDirection()[2], trackID, parentID, isIncident, isExit, track->GetVertexKineticEnergy(), edep, niel, pdgid, charge, creatorProcess); } } } if(pvname == "ddEpi2" && histoManager) { if(creatorProcess!=-1) { histoManager->FillTree("radioactiveDecay_ddEpi2", rndmSeed, evtNumber, x,y,z,momentum.x(), momentum.y(), momentum.z(), trackID, parentID, isIncident, isExit, kineticEnergy, edep, niel, pdgid, charge, creatorProcess); G4cout << "evt number = " << evtNumber << G4endl; G4cout << track->GetPosition() << G4endl; G4cout << track->GetVertexPosition() << G4endl << G4endl; if(radioactiveTrackIDEpi2!=trackID) // only want to count once the source of the decay { radioactiveTrackIDEpi2 = trackID; histoManager->FillTree("radioactiveDecay_ddEpi2_Vertex", rndmSeed, evtNumber, track->GetVertexPosition()[0],track->GetVertexPosition()[1],track->GetVertexPosition()[2], track->GetVertexMomentumDirection()[0],track->GetVertexMomentumDirection()[1],track->GetVertexMomentumDirection()[2], trackID, parentID, isIncident, isExit, track->GetVertexKineticEnergy(), edep, niel, pdgid, charge, creatorProcess); } } } if(pvname == "ddSub1" && pdgid==2212) { G4double e = step->GetTrack()->GetKineticEnergy(); G4double de = step->GetTotalEnergyDeposit()/keV; G4double dx = step->GetStepLength()/um; /// G4cout << e << " " << de << " " << dx << G4endl; if(e>5 && e<60 && de>0 && dx>0) histoManager->FillProfile("letscan",e, de/dx); } */ if(pvname == "Epi" || pvname == "Substrate") RangeTelescopeSteppingAction(step); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4double SteppingAction::GetThickness(const G4Step* step) { G4Box* box = (G4Box*)step->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetSolid(); G4double t = box->GetZHalfLength()*2; //delete box; return t/um; } void SteppingAction::TrackerSteppingAction(const G4Step* step) { //Can be used to write step by step particle track info to output ascii file } void SteppingAction::RangeTelescopeSteppingAction(const G4Step* step) { // G4cout << "RangeTelescopeSteppingAction::Test" << G4endl; //G4double EpiThick=14; // G4double SubThick=735; //G4double lambdaSub=4.68; G4StepPoint* preStepPoint = step->GetPreStepPoint(); G4String vol = preStepPoint->GetPhysicalVolume()->GetName(); G4ThreeVector momentum= preStepPoint->GetMomentumDirection(); G4TouchableHistory*theTouchable = (G4TouchableHistory*)(preStepPoint->GetTouchable()); G4String name = step->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName(); G4int PartNum = step->GetTrack()->GetDynamicParticle()->GetDefinition()->GetPDGEncoding(); G4double e = step ->GetPostStepPoint()->GetKineticEnergy(); G4String process = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); G4double edep = step->GetTotalEnergyDeposit(); G4double niel = step->GetNonIonizingEnergyDeposit(); G4ThreeVector worldPosition = preStepPoint->GetPosition(); G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition); G4double x = localPosition.x(); G4double y = localPosition.y(); G4double z = localPosition.z(); ///G4double xAbs= step->GetPostStepPoint()->GetPosition().x(); G4int trackID= step->GetTrack()->GetTrackID(); G4int parentID = step->GetTrack()->GetParentID(); G4int motherCopyNo = theTouchable->GetReplicaNumber(1); // get depth==1 to get mother replica number in case ///motherCopyNoStrip = tracker->FindStripCopyNumber(xAbs); // G4cout << "RangeTelescopeSteppingAction::Test trackID = " << trackID << G4endl; // G4cout << "RangeTelescopeSteppingAction::Test pdgid = " << PartNum << G4endl; // G4cout << "RangeTelescopeSteppingAction::Test z = " << worldPosition.z() << G4endl; // G4cout << "RangeTelescopeSteppingAction::Test replica number = " << theTouchable->GetReplicaNumber(1) << G4endl; // G4cout << "RangeTelescopeSteppingAction::Test process = " << process << G4endl; //HistoManager* histoManager = HistoManager::GetHistoManager(); if(histoManager) { if(trackID==1) histoManager->FillHisto2D("primaryCopyVsZ", worldPosition.z(), motherCopyNo); else histoManager->FillHisto2D("secondaryCopyVsZ", worldPosition.z(), motherCopyNo); } // truith plane of the RT which will take all incident particles if(motherCopyNo==0 && preStepPoint->GetStepStatus() == fGeomBoundary && vol=="Epi") { G4int evtNumber = 0; const G4Event* evt = G4RunManager::GetRunManager()->GetCurrentEvent(); if(evt) evtNumber = evt->GetEventID(); G4double rndmSeed = G4Random::getTheSeed(); if(histoManager) { histoManager->FillTree("RT_TruthPlane", rndmSeed, evtNumber, worldPosition.x(), worldPosition.y(), worldPosition.z(), momentum.x(), momentum.y(), momentum.z(), trackID, parentID, 1, 0, e, edep, niel, PartNum, -99, -99); } fEventAction->SetRTTruthEnergy(e); } //Energy deposited by ionization in epitaxial and substrate is passed to trasnport classes for charge sharing // pass the motherCopyNo to EventAction to store last layer per event if(vol=="Epi" && PartNum==2212) { fEventAction->SetLastRTLayerPerEvent(motherCopyNo); } if (process=="eIoni"|| process=="hIoni") { //G4cout << "SteppingAction edep before = " << edep << G4endl; if(vol=="Epi") { if(EpiThick < 0) EpiThick = GetThickness(step); //G4cout << "EpiThick = " << EpiThick << G4endl; fEventAction->AddEpiEdep(motherCopyNo, edep); if (z/um+(EpiThick/2)<1.1) { //Lower collection efficiency when diode is placed edep = edep*0.8; if( std::isfinite(edep) ) //check for NaN and inf fTrackingAction -> GetTransport()-> Gauss(edep,x,y,z,motherCopyNo,"CMOS"); } else { if( std::isfinite(edep) ) //check for NaN and inf fTrackingAction -> GetTransport()-> Gauss(edep,x,y,z,motherCopyNo,"CMOS"); } } if(vol=="Substrate") { if(SubThick < 0) SubThick = GetThickness(step); //G4cout << "SubThick = " << SubThick << G4endl; //Charge collection in substrate G4double lambdaSub=4.68; G4double charge= exp(-(SubThick-z/um)/lambdaSub); if(charge>0.01) //if(std::isfinite(charge) && charge>0.01 && charge<1 ) { //just not carrying on amount of charge which is too small edep = edep * charge; if( std::isfinite(edep) ) //check for NaN and inf //G4cout << "SteppingAction charge sub = " << charge << G4endl; fTrackingAction -> GetTransport()-> Gauss(edep,x,y,z,motherCopyNo,"CMOS"); } } } } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // note values from paper in fGy/m2 Double_t GetNeutronKermaCoefficient(Double_t energy) { double coeff = 0.0; if(energy < 0.13) coeff = 0.00161; else if(energy < 0.14) coeff = 0.00110; else if(energy < 0.15) coeff = 0.000728; else if(energy < 0.16) coeff = 0.00222; else if(energy < 0.17) coeff = 0.00984; else if(energy < 0.18) coeff = 0.0335; else if(energy < 0.19) coeff = 0.0595; else if(energy < 0.22) coeff = 0.0449; else if(energy < 0.24) coeff = 0.0356; else if(energy < 0.26) coeff = 0.0310; else if(energy < 0.28) coeff = 0.0289; else if(energy < 0.30) coeff = 0.0280; else if(energy < 0.32) coeff = 0.0276; else if(energy < 0.34) coeff = 0.0273; else if(energy < 0.36) coeff = 0.0279; else if(energy < 0.38) coeff = 0.0279; else if(energy < 0.40) coeff = 0.0295; else if(energy < 0.44) coeff = 0.0298; else if(energy < 0.48) coeff = 0.0311; else if(energy < 0.52) coeff = 0.0332; else if(energy < 0.56) coeff = 0.0500; else if(energy < 0.60) coeff = 0.0514; else if(energy < 0.88) coeff = 0.0484; else if(energy < 1.05) coeff = 0.0480; else if(energy < 1.40) coeff = 0.0623; else if(energy < 1.70) coeff = 0.127; else if(energy < 1.90) coeff = 0.0950; else if(energy < 2.40) coeff = 0.0812; else if(energy < 2.60) coeff = 0.100; else if(energy < 3.00) coeff = 0.100; else if(energy < 3.20) coeff = 0.112; else if(energy < 3.60) coeff = 0.103; else if(energy < 3.80) coeff = 0.0862; else if(energy < 5.00) coeff = 0.175; else if(energy < 6.00) coeff = 0.218; else if(energy < 7.00) coeff = 0.413; else if(energy < 8.00) coeff = 0.719; else if(energy < 9.00) coeff = 0.769; else if(energy < 10.00) coeff = 0.949; else if(energy < 12.00) coeff = 1.05; else if(energy < 13.00) coeff = 1.14; else if(energy < 14.00) coeff = 1.21; else if(energy < 15.00) coeff = 1.28; else if(energy < 16.00) coeff = 1.33; else if(energy < 18.00) coeff = 1.39; else if(energy < 20.00) coeff = 1.45; else if(energy < 23.00) coeff = 1.61; else if(energy < 25.00) coeff = 1.74; else if(energy < 27.00) coeff = 1.87; else if(energy < 29.00) coeff = 2.00; else if(energy < 31.00) coeff = 2.14; else coeff = 2.43; //convert from fGy/m2 to fGy/cm2 coeff = coeff * 1e4; //convert from fGy/cm2 to Gy/cm2 coeff = coeff * 1e-15; //return value in Gy/neutron/cm2 return coeff; } G4double SteppingAction::GetProton1MeVNeutronEquivalence(G4double energy) { // for now use my power law and update when get proper data //double pow = -0.52; //double a = 14.00; double neq(0);// = a * TMath::Power(energy, pow); //std::cout << energy << " " << neq << std::endl; if(energy< 1.000E-03 ) neq = 7.590E+03 ; else if(energy< 2.000E-03 ) neq = 5.170E+03 ; else if(energy< 3.000E-03 ) neq = 3.976E+03 ; else if(energy< 5.000E-03 ) neq = 2.778E+03 ; else if(energy< 7.000E-03 ) neq = 2.167E+03 ; else if(energy< 1.000E-02 ) neq = 1.650E+03 ; else if(energy< 2.000E-02 ) neq = 9.515E+02 ; else if(energy< 3.000E-02 ) neq = 6.824E+02 ; else if(energy< 5.000E-02 ) neq = 4.447E+02 ; else if(energy< 7.000E-02 ) neq = 3.338E+02 ; else if(energy< 1.000E-01 ) neq = 2.456E+02 ; else if(energy< 2.000E-01 ) neq = 1.336E+02 ; else if(energy< 3.000E-01 ) neq = 9.127E+01 ; else if(energy< 5.000E-01 ) neq = 5.862E+01 ; else if(energy< 7.000E-01 ) neq = 4.307E+01 ; else if(energy< 1.000E+00 ) neq = 3.133E+01 ; else if(energy< 2.000E+00 ) neq = 1.618E+01 ; else if(energy< 3.000E+00 ) neq = 1.101E+01 ; else if(energy< 5.000E+00 ) neq = 6.756E+00 ; else if(energy< 7.000E+00 ) neq = 5.140E+00 ; else if(energy< 1.000E+01 ) neq = 3.871E+00 ; else if(energy< 2.000E+01 ) neq = 2.632E+00 ; else if(energy< 3.000E+01 ) neq = 2.346E+00 ; else if(energy< 5.000E+01 ) neq = 1.907E+00 ; else if(energy< 7.000E+01 ) neq = 1.552E+00 ; else if(energy< 1.000E+02 ) neq = 1.276E+00 ; else if(energy< 1.050E+02 ) neq = 1.244E+00 ; else if(energy< 1.150E+02 ) neq = 1.202E+00 ; else if(energy< 1.250E+02 ) neq = 1.161E+00 ; else if(energy< 1.350E+02 ) neq = 1.124E+00 ; else if(energy< 1.450E+02 ) neq = 1.099E+00 ; else if(energy< 1.550E+02 ) neq = 1.075E+00 ; else if(energy< 1.650E+02 ) neq = 1.051E+00 ; else if(energy< 1.750E+02 ) neq = 1.028E+00 ; else if(energy< 1.850E+02 ) neq = 1.013E+00 ; else if(energy< 1.950E+02 ) neq = 9.968E-01 ; else if(energy< 2.050E+02 ) neq = 9.811E-01 ; return neq; } G4double SteppingAction::GetNeutron1MeVNeutronEquivalence(G4double energy) { double neq(0); if(energy< 1.025E-10 ) neq = 1.575E-02 ; else if(energy < 1.075E-10 ) neq = 1.537E-02 ; else if(energy < 1.125E-10 ) neq = 1.503E-02 ; else if(energy < 1.175E-10 ) neq = 1.470E-02 ; else if(energy < 1.238E-10 ) neq = 1.432E-02 ; else if(energy < 1.313E-10 ) neq = 1.391E-02 ; else if(energy < 1.388E-10 ) neq = 1.353E-02 ; else if(energy < 1.463E-10 ) neq = 1.317E-02 ; else if(energy < 1.550E-10 ) neq = 1.280E-02 ; else if(energy < 1.650E-10 ) neq = 1.242E-02 ; else if(energy < 1.750E-10 ) neq = 1.206E-02 ; else if(energy < 1.850E-10 ) neq = 1.172E-02 ; else if(energy < 1.950E-10 ) neq = 1.142E-02 ; else if(energy < 2.050E-10 ) neq = 1.113E-02 ; else if(energy < 2.150E-10 ) neq = 1.087E-02 ; else if(energy < 2.250E-10 ) neq = 1.063E-02 ; else if(energy < 2.350E-10 ) neq = 1.039E-02 ; else if(energy < 2.475E-10 ) neq = 1.013E-02 ; else if(energy < 2.625E-10 ) neq = 9.834E-03 ; else if(energy < 2.750E-10 ) neq = 9.609E-03 ; else if(energy < 2.900E-10 ) neq = 9.356E-03 ; else if(energy < 3.100E-10 ) neq = 9.048E-03 ; else if(energy < 3.300E-10 ) neq = 8.770E-03 ; else if(energy < 3.500E-10 ) neq = 8.516E-03 ; else if(energy < 3.700E-10 ) neq = 8.290E-03 ; else if(energy < 3.900E-10 ) neq = 8.075E-03 ; else if(energy < 4.125E-10 ) neq = 7.847E-03 ; else if(energy < 4.375E-10 ) neq = 7.622E-03 ; else if(energy < 4.625E-10 ) neq = 7.410E-03 ; else if(energy < 4.875E-10 ) neq = 7.219E-03 ; else if(energy < 5.125E-10 ) neq = 7.039E-03 ; else if(energy < 5.375E-10 ) neq = 6.875E-03 ; else if(energy < 5.625E-10 ) neq = 6.718E-03 ; else if(energy < 5.875E-10 ) neq = 6.575E-03 ; else if(energy < 6.150E-10 ) neq = 6.425E-03 ; else if(energy < 6.450E-10 ) neq = 6.274E-03 ; else if(energy < 6.750E-10 ) neq = 6.134E-03 ; else if(energy < 7.050E-10 ) neq = 6.002E-03 ; else if(energy < 7.400E-10 ) neq = 5.864E-03 ; else if(energy < 7.800E-10 ) neq = 5.709E-03 ; else if(energy < 8.200E-10 ) neq = 5.567E-03 ; else if(energy < 8.600E-10 ) neq = 5.438E-03 ; else if(energy < 9.000E-10 ) neq = 5.312E-03 ; else if(energy < 9.400E-10 ) neq = 5.200E-03 ; else if(energy < 9.800E-10 ) neq = 5.091E-03 ; else if(energy < 1.025E-09 ) neq = 4.978E-03 ; else if(energy < 1.075E-09 ) neq = 4.861E-03 ; else if(energy < 1.125E-09 ) neq = 4.751E-03 ; else if(energy < 1.175E-09 ) neq = 4.649E-03 ; else if(energy < 1.238E-09 ) neq = 4.530E-03 ; else if(energy < 1.313E-09 ) neq = 4.398E-03 ; else if(energy < 1.388E-09 ) neq = 4.278E-03 ; else if(energy < 1.463E-09 ) neq = 4.171E-03 ; else if(energy < 1.550E-09 ) neq = 4.049E-03 ; else if(energy < 1.650E-09 ) neq = 3.925E-03 ; else if(energy < 1.750E-09 ) neq = 3.810E-03 ; else if(energy < 1.850E-09 ) neq = 3.706E-03 ; else if(energy < 1.950E-09 ) neq = 3.609E-03 ; else if(energy < 2.050E-09 ) neq = 3.520E-03 ; else if(energy < 2.150E-09 ) neq = 3.437E-03 ; else if(energy < 2.250E-09 ) neq = 3.360E-03 ; else if(energy < 2.350E-09 ) neq = 3.287E-03 ; else if(energy < 2.475E-09 ) neq = 3.203E-03 ; else if(energy < 2.625E-09 ) neq = 3.110E-03 ; else if(energy < 2.750E-09 ) neq = 3.039E-03 ; else if(energy < 2.900E-09 ) neq = 2.962E-03 ; else if(energy < 3.100E-09 ) neq = 2.864E-03 ; else if(energy < 3.300E-09 ) neq = 2.776E-03 ; else if(energy < 3.500E-09 ) neq = 2.695E-03 ; else if(energy < 3.700E-09 ) neq = 2.621E-03 ; else if(energy < 3.900E-09 ) neq = 2.552E-03 ; else if(energy < 4.125E-09 ) neq = 2.482E-03 ; else if(energy < 4.375E-09 ) neq = 2.409E-03 ; else if(energy < 4.625E-09 ) neq = 2.344E-03 ; else if(energy < 4.875E-09 ) neq = 2.283E-03 ; else if(energy < 5.125E-09 ) neq = 2.226E-03 ; else if(energy < 5.375E-09 ) neq = 2.174E-03 ; else if(energy < 5.625E-09 ) neq = 2.125E-03 ; else if(energy < 5.875E-09 ) neq = 2.082E-03 ; else if(energy < 6.150E-09 ) neq = 2.034E-03 ; else if(energy < 6.450E-09 ) neq = 1.985E-03 ; else if(energy < 6.750E-09 ) neq = 1.941E-03 ; else if(energy < 7.050E-09 ) neq = 1.898E-03 ; else if(energy < 7.400E-09 ) neq = 1.853E-03 ; else if(energy < 7.800E-09 ) neq = 1.805E-03 ; else if(energy < 8.200E-09 ) neq = 1.760E-03 ; else if(energy < 8.600E-09 ) neq = 1.718E-03 ; else if(energy < 9.000E-09 ) neq = 1.680E-03 ; else if(energy < 9.400E-09 ) neq = 1.644E-03 ; else if(energy < 9.800E-09 ) neq = 1.610E-03 ; else if(energy < 1.025E-08 ) neq = 1.574E-03 ; else if(energy < 1.075E-08 ) neq = 1.537E-03 ; else if(energy < 1.125E-08 ) neq = 1.503E-03 ; else if(energy < 1.175E-08 ) neq = 1.472E-03 ; else if(energy < 1.238E-08 ) neq = 1.433E-03 ; else if(energy < 1.313E-08 ) neq = 1.392E-03 ; else if(energy < 1.388E-08 ) neq = 1.353E-03 ; else if(energy < 1.463E-08 ) neq = 1.318E-03 ; else if(energy < 1.550E-08 ) neq = 1.280E-03 ; else if(energy < 1.650E-08 ) neq = 1.241E-03 ; else if(energy < 1.750E-08 ) neq = 1.204E-03 ; else if(energy < 1.850E-08 ) neq = 1.172E-03 ; else if(energy < 1.950E-08 ) neq = 1.141E-03 ; else if(energy < 2.050E-08 ) neq = 1.113E-03 ; else if(energy < 2.150E-08 ) neq = 1.087E-03 ; else if(energy < 2.250E-08 ) neq = 1.062E-03 ; else if(energy < 2.350E-08 ) neq = 1.040E-03 ; else if(energy < 2.475E-08 ) neq = 1.013E-03 ; else if(energy < 2.625E-08 ) neq = 9.841E-04 ; else if(energy < 2.750E-08 ) neq = 9.619E-04 ; else if(energy < 2.900E-08 ) neq = 9.359E-04 ; else if(energy < 3.100E-08 ) neq = 9.058E-04 ; else if(energy < 3.300E-08 ) neq = 8.773E-04 ; else if(energy < 3.500E-08 ) neq = 8.523E-04 ; else if(energy < 3.700E-08 ) neq = 8.285E-04 ; else if(energy < 3.900E-08 ) neq = 8.073E-04 ; else if(energy < 4.125E-08 ) neq = 7.847E-04 ; else if(energy < 4.375E-08 ) neq = 7.619E-04 ; else if(energy < 4.625E-08 ) neq = 7.410E-04 ; else if(energy < 4.875E-08 ) neq = 7.215E-04 ; else if(energy < 5.125E-08 ) neq = 7.036E-04 ; else if(energy < 5.375E-08 ) neq = 6.871E-04 ; else if(energy < 5.625E-08 ) neq = 6.716E-04 ; else if(energy < 5.875E-08 ) neq = 6.578E-04 ; else if(energy < 6.150E-08 ) neq = 6.428E-04 ; else if(energy < 6.450E-08 ) neq = 6.273E-04 ; else if(energy < 6.750E-08 ) neq = 6.136E-04 ; else if(energy < 7.050E-08 ) neq = 6.001E-04 ; else if(energy < 7.400E-08 ) neq = 5.861E-04 ; else if(energy < 7.800E-08 ) neq = 5.707E-04 ; else if(energy < 8.200E-08 ) neq = 5.567E-04 ; else if(energy < 8.600E-08 ) neq = 5.434E-04 ; else if(energy < 9.000E-08 ) neq = 5.313E-04 ; else if(energy < 9.400E-08 ) neq = 5.197E-04 ; else if(energy < 9.800E-08 ) neq = 5.090E-04 ; else if(energy < 1.025E-07 ) neq = 4.976E-04 ; else if(energy < 1.075E-07 ) neq = 4.864E-04 ; else if(energy < 1.125E-07 ) neq = 4.753E-04 ; else if(energy < 1.175E-07 ) neq = 4.646E-04 ; else if(energy < 1.238E-07 ) neq = 4.531E-04 ; else if(energy < 1.313E-07 ) neq = 4.397E-04 ; else if(energy < 1.388E-07 ) neq = 4.279E-04 ; else if(energy < 1.463E-07 ) neq = 4.166E-04 ; else if(energy < 1.550E-07 ) neq = 4.049E-04 ; else if(energy < 1.650E-07 ) neq = 3.924E-04 ; else if(energy < 1.750E-07 ) neq = 3.810E-04 ; else if(energy < 1.850E-07 ) neq = 3.705E-04 ; else if(energy < 1.950E-07 ) neq = 3.607E-04 ; else if(energy < 2.050E-07 ) neq = 3.518E-04 ; else if(energy < 2.150E-07 ) neq = 3.437E-04 ; else if(energy < 2.250E-07 ) neq = 3.362E-04 ; else if(energy < 2.350E-07 ) neq = 3.286E-04 ; else if(energy < 2.475E-07 ) neq = 3.204E-04 ; else if(energy < 2.625E-07 ) neq = 3.111E-04 ; else if(energy < 2.750E-07 ) neq = 3.039E-04 ; else if(energy < 2.900E-07 ) neq = 2.962E-04 ; else if(energy < 3.100E-07 ) neq = 2.865E-04 ; else if(energy < 3.300E-07 ) neq = 2.775E-04 ; else if(energy < 3.500E-07 ) neq = 2.695E-04 ; else if(energy < 3.700E-07 ) neq = 2.620E-04 ; else if(energy < 3.900E-07 ) neq = 2.551E-04 ; else if(energy < 4.125E-07 ) neq = 2.483E-04 ; else if(energy < 4.375E-07 ) neq = 2.410E-04 ; else if(energy < 4.625E-07 ) neq = 2.343E-04 ; else if(energy < 4.875E-07 ) neq = 2.284E-04 ; else if(energy < 5.125E-07 ) neq = 2.226E-04 ; else if(energy < 5.375E-07 ) neq = 2.176E-04 ; else if(energy < 5.625E-07 ) neq = 2.128E-04 ; else if(energy < 5.875E-07 ) neq = 2.082E-04 ; else if(energy < 6.150E-07 ) neq = 2.034E-04 ; else if(energy < 6.450E-07 ) neq = 1.985E-04 ; else if(energy < 6.750E-07 ) neq = 1.941E-04 ; else if(energy < 7.050E-07 ) neq = 1.898E-04 ; else if(energy < 7.400E-07 ) neq = 1.852E-04 ; else if(energy < 7.800E-07 ) neq = 1.805E-04 ; else if(energy < 8.200E-07 ) neq = 1.762E-04 ; else if(energy < 8.600E-07 ) neq = 1.718E-04 ; else if(energy < 9.000E-07 ) neq = 1.679E-04 ; else if(energy < 9.400E-07 ) neq = 1.645E-04 ; else if(energy < 9.800E-07 ) neq = 1.610E-04 ; else if(energy < 1.025E-06 ) neq = 1.575E-04 ; else if(energy < 1.075E-06 ) neq = 1.542E-04 ; else if(energy < 1.125E-06 ) neq = 1.508E-04 ; else if(energy < 1.175E-06 ) neq = 1.474E-04 ; else if(energy < 1.238E-06 ) neq = 1.434E-04 ; else if(energy < 1.313E-06 ) neq = 1.393E-04 ; else if(energy < 1.388E-06 ) neq = 1.356E-04 ; else if(energy < 1.463E-06 ) neq = 1.318E-04 ; else if(energy < 1.550E-06 ) neq = 1.281E-04 ; else if(energy < 1.650E-06 ) neq = 1.242E-04 ; else if(energy < 1.750E-06 ) neq = 1.204E-04 ; else if(energy < 1.850E-06 ) neq = 1.172E-04 ; else if(energy < 1.950E-06 ) neq = 1.141E-04 ; else if(energy < 2.050E-06 ) neq = 1.114E-04 ; else if(energy < 2.150E-06 ) neq = 1.090E-04 ; else if(energy < 2.250E-06 ) neq = 1.066E-04 ; else if(energy < 2.350E-06 ) neq = 1.042E-04 ; else if(energy < 2.475E-06 ) neq = 1.013E-04 ; else if(energy < 2.625E-06 ) neq = 9.853E-05 ; else if(energy < 2.750E-06 ) neq = 9.627E-05 ; else if(energy < 2.900E-06 ) neq = 9.360E-05 ; else if(energy < 3.100E-06 ) neq = 9.056E-05 ; else if(energy < 3.300E-06 ) neq = 8.776E-05 ; else if(energy < 3.500E-06 ) neq = 8.513E-05 ; else if(energy < 3.700E-06 ) neq = 8.287E-05 ; else if(energy < 3.900E-06 ) neq = 8.065E-05 ; else if(energy < 4.125E-06 ) neq = 7.860E-05 ; else if(energy < 4.375E-06 ) neq = 7.644E-05 ; else if(energy < 4.625E-06 ) neq = 7.428E-05 ; else if(energy < 4.875E-06 ) neq = 7.217E-05 ; else if(energy < 5.125E-06 ) neq = 7.045E-05 ; else if(energy < 5.375E-06 ) neq = 6.885E-05 ; else if(energy < 5.625E-06 ) neq = 6.725E-05 ; else if(energy < 5.875E-06 ) neq = 6.569E-05 ; else if(energy < 6.150E-06 ) neq = 6.427E-05 ; else if(energy < 6.450E-06 ) neq = 6.278E-05 ; else if(energy < 6.750E-06 ) neq = 6.129E-05 ; else if(energy < 7.050E-06 ) neq = 5.998E-05 ; else if(energy < 7.400E-06 ) neq = 5.858E-05 ; else if(energy < 7.800E-06 ) neq = 5.701E-05 ; else if(energy < 8.200E-06 ) neq = 5.572E-05 ; else if(energy < 8.600E-06 ) neq = 5.450E-05 ; else if(energy < 9.000E-06 ) neq = 5.327E-05 ; else if(energy < 9.400E-06 ) neq = 5.205E-05 ; else if(energy < 9.800E-06 ) neq = 5.087E-05 ; else if(energy < 1.025E-05 ) neq = 4.981E-05 ; else if(energy < 1.075E-05 ) neq = 4.868E-05 ; else if(energy < 1.125E-05 ) neq = 4.754E-05 ; else if(energy < 1.175E-05 ) neq = 4.644E-05 ; else if(energy < 1.238E-05 ) neq = 4.531E-05 ; else if(energy < 1.313E-05 ) neq = 4.399E-05 ; else if(energy < 1.388E-05 ) neq = 4.274E-05 ; else if(energy < 1.463E-05 ) neq = 4.167E-05 ; else if(energy < 1.550E-05 ) neq = 4.045E-05 ; else if(energy < 1.650E-05 ) neq = 3.929E-05 ; else if(energy < 1.750E-05 ) neq = 3.821E-05 ; else if(energy < 1.850E-05 ) neq = 3.712E-05 ; else if(energy < 1.950E-05 ) neq = 3.607E-05 ; else if(energy < 2.050E-05 ) neq = 3.522E-05 ; else if(energy < 2.150E-05 ) neq = 3.442E-05 ; else if(energy < 2.250E-05 ) neq = 3.361E-05 ; else if(energy < 2.350E-05 ) neq = 3.284E-05 ; else if(energy < 2.475E-05 ) neq = 3.204E-05 ; else if(energy < 2.625E-05 ) neq = 3.110E-05 ; else if(energy < 2.750E-05 ) neq = 3.035E-05 ; else if(energy < 2.900E-05 ) neq = 2.959E-05 ; else if(energy < 3.100E-05 ) neq = 2.860E-05 ; else if(energy < 3.300E-05 ) neq = 2.778E-05 ; else if(energy < 3.500E-05 ) neq = 2.702E-05 ; else if(energy < 3.700E-05 ) neq = 2.625E-05 ; else if(energy < 3.900E-05 ) neq = 2.550E-05 ; else if(energy < 4.125E-05 ) neq = 2.483E-05 ; else if(energy < 4.375E-05 ) neq = 2.412E-05 ; else if(energy < 4.625E-05 ) neq = 2.342E-05 ; else if(energy < 4.875E-05 ) neq = 2.282E-05 ; else if(energy < 5.125E-05 ) neq = 2.227E-05 ; else if(energy < 5.375E-05 ) neq = 2.172E-05 ; else if(energy < 5.625E-05 ) neq = 2.123E-05 ; else if(energy < 5.875E-05 ) neq = 2.079E-05 ; else if(energy < 6.150E-05 ) neq = 2.030E-05 ; else if(energy < 6.450E-05 ) neq = 1.985E-05 ; else if(energy < 6.750E-05 ) neq = 1.944E-05 ; else if(energy < 7.050E-05 ) neq = 1.904E-05 ; else if(energy < 7.400E-05 ) neq = 1.856E-05 ; else if(energy < 7.800E-05 ) neq = 1.803E-05 ; else if(energy < 8.200E-05 ) neq = 1.761E-05 ; else if(energy < 8.600E-05 ) neq = 1.721E-05 ; else if(energy < 9.000E-05 ) neq = 1.681E-05 ; else if(energy < 9.400E-05 ) neq = 1.642E-05 ; else if(energy < 9.800E-05 ) neq = 1.610E-05 ; else if(energy < 1.025E-04 ) neq = 1.575E-05 ; else if(energy < 1.075E-04 ) neq = 1.536E-05 ; else if(energy < 1.125E-04 ) neq = 1.502E-05 ; else if(energy < 1.175E-04 ) neq = 1.470E-05 ; else if(energy < 1.238E-04 ) neq = 1.432E-05 ; else if(energy < 1.313E-04 ) neq = 1.393E-05 ; else if(energy < 1.388E-04 ) neq = 1.357E-05 ; else if(energy < 1.463E-04 ) neq = 1.321E-05 ; else if(energy < 1.550E-04 ) neq = 1.714E-05 ; else if(energy < 1.650E-04 ) neq = 6.792E-05 ; else if(energy < 1.750E-04 ) neq = 1.315E-04 ; else if(energy < 1.850E-04 ) neq = 1.938E-04 ; else if(energy < 1.950E-04 ) neq = 2.250E-04 ; else if(energy < 2.050E-04 ) neq = 2.433E-04 ; else if(energy < 2.150E-04 ) neq = 2.614E-04 ; else if(energy < 2.250E-04 ) neq = 2.754E-04 ; else if(energy < 2.350E-04 ) neq = 2.867E-04 ; else if(energy < 2.475E-04 ) neq = 3.007E-04 ; else if(energy < 2.625E-04 ) neq = 3.177E-04 ; else if(energy < 2.750E-04 ) neq = 3.318E-04 ; else if(energy < 2.900E-04 ) neq = 3.486E-04 ; else if(energy < 3.100E-04 ) neq = 3.712E-04 ; else if(energy < 3.300E-04 ) neq = 3.937E-04 ; else if(energy < 3.500E-04 ) neq = 4.163E-04 ; else if(energy < 3.700E-04 ) neq = 4.388E-04 ; else if(energy < 3.900E-04 ) neq = 4.613E-04 ; else if(energy < 4.125E-04 ) neq = 4.865E-04 ; else if(energy < 4.375E-04 ) neq = 5.146E-04 ; else if(energy < 4.625E-04 ) neq = 5.426E-04 ; else if(energy < 4.875E-04 ) neq = 5.706E-04 ; else if(energy < 5.125E-04 ) neq = 5.985E-04 ; else if(energy < 5.375E-04 ) neq = 6.264E-04 ; else if(energy < 5.625E-04 ) neq = 6.543E-04 ; else if(energy < 5.875E-04 ) neq = 6.822E-04 ; else if(energy < 6.150E-04 ) neq = 7.128E-04 ; else if(energy < 6.450E-04 ) neq = 7.462E-04 ; else if(energy < 6.750E-04 ) neq = 7.795E-04 ; else if(energy < 7.050E-04 ) neq = 8.128E-04 ; else if(energy < 7.400E-04 ) neq = 8.515E-04 ; else if(energy < 7.800E-04 ) neq = 8.957E-04 ; else if(energy < 8.200E-04 ) neq = 9.399E-04 ; else if(energy < 8.600E-04 ) neq = 9.840E-04 ; else if(energy < 9.000E-04 ) neq = 1.028E-03 ; else if(energy < 9.400E-04 ) neq = 1.072E-03 ; else if(energy < 9.800E-04 ) neq = 1.116E-03 ; else if(energy < 1.025E-03 ) neq = 1.165E-03 ; else if(energy < 1.075E-03 ) neq = 1.220E-03 ; else if(energy < 1.125E-03 ) neq = 1.274E-03 ; else if(energy < 1.175E-03 ) neq = 1.329E-03 ; else if(energy < 1.238E-03 ) neq = 1.397E-03 ; else if(energy < 1.313E-03 ) neq = 1.478E-03 ; else if(energy < 1.388E-03 ) neq = 1.559E-03 ; else if(energy < 1.463E-03 ) neq = 1.640E-03 ; else if(energy < 1.550E-03 ) neq = 1.734E-03 ; else if(energy < 1.650E-03 ) neq = 1.842E-03 ; else if(energy < 1.750E-03 ) neq = 1.949E-03 ; else if(energy < 1.850E-03 ) neq = 2.056E-03 ; else if(energy < 1.950E-03 ) neq = 2.162E-03 ; else if(energy < 2.050E-03 ) neq = 2.269E-03 ; else if(energy < 2.150E-03 ) neq = 2.382E-03 ; else if(energy < 2.250E-03 ) neq = 4.405E-03 ; else if(energy < 2.350E-03 ) neq = 2.589E-03 ; else if(energy < 2.475E-03 ) neq = 2.717E-03 ; else if(energy < 2.625E-03 ) neq = 2.874E-03 ; else if(energy < 2.750E-03 ) neq = 3.004E-03 ; else if(energy < 2.900E-03 ) neq = 3.159E-03 ; else if(energy < 3.100E-03 ) neq = 3.366E-03 ; else if(energy < 3.300E-03 ) neq = 3.573E-03 ; else if(energy < 3.500E-03 ) neq = 3.778E-03 ; else if(energy < 3.700E-03 ) neq = 3.982E-03 ; else if(energy < 3.900E-03 ) neq = 4.186E-03 ; else if(energy < 4.125E-03 ) neq = 4.413E-03 ; else if(energy < 4.375E-03 ) neq = 4.666E-03 ; else if(energy < 4.625E-03 ) neq = 4.917E-03 ; else if(energy < 4.875E-03 ) neq = 5.492E-03 ; else if(energy < 5.125E-03 ) neq = 5.419E-03 ; else if(energy < 5.375E-03 ) neq = 5.662E-03 ; else if(energy < 5.625E-03 ) neq = 5.907E-03 ; else if(energy < 5.875E-03 ) neq = 6.151E-03 ; else if(energy < 6.150E-03 ) neq = 6.418E-03 ; else if(energy < 6.450E-03 ) neq = 6.710E-03 ; else if(energy < 6.750E-03 ) neq = 6.999E-03 ; else if(energy < 7.050E-03 ) neq = 7.285E-03 ; else if(energy < 7.400E-03 ) neq = 7.618E-03 ; else if(energy < 7.800E-03 ) neq = 7.998E-03 ; else if(energy < 8.200E-03 ) neq = 8.372E-03 ; else if(energy < 8.600E-03 ) neq = 8.745E-03 ; else if(energy < 9.000E-03 ) neq = 9.115E-03 ; else if(energy < 9.400E-03 ) neq = 9.477E-03 ; else if(energy < 9.800E-03 ) neq = 9.838E-03 ; else if(energy < 1.025E-02 ) neq = 1.024E-02 ; else if(energy < 1.075E-02 ) neq = 1.069E-02 ; else if(energy < 1.125E-02 ) neq = 1.115E-02 ; else if(energy < 1.175E-02 ) neq = 1.159E-02 ; else if(energy < 1.238E-02 ) neq = 1.213E-02 ; else if(energy < 1.313E-02 ) neq = 1.277E-02 ; else if(energy < 1.388E-02 ) neq = 1.340E-02 ; else if(energy < 1.463E-02 ) neq = 1.403E-02 ; else if(energy < 1.550E-02 ) neq = 1.557E-02 ; else if(energy < 1.650E-02 ) neq = 1.554E-02 ; else if(energy < 1.750E-02 ) neq = 1.634E-02 ; else if(energy < 1.850E-02 ) neq = 1.710E-02 ; else if(energy < 1.950E-02 ) neq = 1.782E-02 ; else if(energy < 2.050E-02 ) neq = 1.854E-02 ; else if(energy < 2.150E-02 ) neq = 1.926E-02 ; else if(energy < 2.250E-02 ) neq = 1.998E-02 ; else if(energy < 2.350E-02 ) neq = 2.069E-02 ; else if(energy < 2.475E-02 ) neq = 2.150E-02 ; else if(energy < 2.625E-02 ) neq = 2.245E-02 ; else if(energy < 2.750E-02 ) neq = 2.317E-02 ; else if(energy < 2.900E-02 ) neq = 2.400E-02 ; else if(energy < 3.100E-02 ) neq = 2.509E-02 ; else if(energy < 3.300E-02 ) neq = 2.599E-02 ; else if(energy < 3.500E-02 ) neq = 2.675E-02 ; else if(energy < 3.700E-02 ) neq = 2.743E-02 ; else if(energy < 3.900E-02 ) neq = 3.247E-02 ; else if(energy < 4.125E-02 ) neq = 2.792E-02 ; else if(energy < 4.375E-02 ) neq = 2.739E-02 ; else if(energy < 4.625E-02 ) neq = 2.598E-02 ; else if(energy < 4.875E-02 ) neq = 2.288E-02 ; else if(energy < 5.125E-02 ) neq = 1.626E-02 ; else if(energy < 5.375E-02 ) neq = 1.485E-02 ; else if(energy < 5.625E-02 ) neq = 5.078E-01 ; else if(energy < 5.875E-02 ) neq = 1.145E-01 ; else if(energy < 6.150E-02 ) neq = 7.793E-02 ; else if(energy < 6.450E-02 ) neq = 6.706E-02 ; else if(energy < 6.750E-02 ) neq = 6.455E-02 ; else if(energy < 7.050E-02 ) neq = 5.988E-02 ; else if(energy < 7.400E-02 ) neq = 5.797E-02 ; else if(energy < 7.800E-02 ) neq = 5.636E-02 ; else if(energy < 8.200E-02 ) neq = 5.494E-02 ; else if(energy < 8.600E-02 ) neq = 5.374E-02 ; else if(energy < 9.000E-02 ) neq = 5.188E-02 ; else if(energy < 9.400E-02 ) neq = 4.992E-02 ; else if(energy < 9.800E-02 ) neq = 4.790E-02 ; else if(energy < 1.025E-01 ) neq = 4.526E-02 ; else if(energy < 1.075E-01 ) neq = 4.162E-02 ; else if(energy < 1.125E-01 ) neq = 3.764E-02 ; else if(energy < 1.175E-01 ) neq = 3.297E-02 ; else if(energy < 1.238E-01 ) neq = 2.645E-02 ; else if(energy < 1.313E-01 ) neq = 1.826E-02 ; else if(energy < 1.388E-01 ) neq = 1.121E-02 ; else if(energy < 1.463E-01 ) neq = 1.090E-02 ; else if(energy < 1.550E-01 ) neq = 4.848E-02 ; else if(energy < 1.650E-01 ) neq = 2.173E-01 ; else if(energy < 1.750E-01 ) neq = 6.987E-01 ; else if(energy < 1.850E-01 ) neq = 1.196E+00 ; else if(energy < 1.950E-01 ) neq = 1.162E+00 ; else if(energy < 2.050E-01 ) neq = 9.710E-01 ; else if(energy < 2.150E-01 ) neq = 8.288E-01 ; else if(energy < 2.250E-01 ) neq = 7.334E-01 ; else if(energy < 2.350E-01 ) neq = 6.709E-01 ; else if(energy < 2.475E-01 ) neq = 6.172E-01 ; else if(energy < 2.625E-01 ) neq = 5.768E-01 ; else if(energy < 2.750E-01 ) neq = 5.547E-01 ; else if(energy < 2.900E-01 ) neq = 5.393E-01 ; else if(energy < 3.100E-01 ) neq = 5.237E-01 ; else if(energy < 3.300E-01 ) neq = 5.198E-01 ; else if(energy < 3.500E-01 ) neq = 5.097E-01 ; else if(energy < 3.700E-01 ) neq = 5.036E-01 ; else if(energy < 3.900E-01 ) neq = 5.320E-01 ; else if(energy < 4.125E-01 ) neq = 5.403E-01 ; else if(energy < 4.375E-01 ) neq = 5.451E-01 ; else if(energy < 4.625E-01 ) neq = 5.581E-01 ; else if(energy < 4.875E-01 ) neq = 5.740E-01 ; else if(energy < 5.125E-01 ) neq = 5.987E-01 ; else if(energy < 5.375E-01 ) neq = 7.554E-01 ; else if(energy < 5.625E-01 ) neq = 1.280E+00 ; else if(energy < 5.875E-01 ) neq = 5.997E-01 ; else if(energy < 6.150E-01 ) neq = 5.425E-01 ; else if(energy < 6.450E-01 ) neq = 5.536E-01 ; else if(energy < 6.750E-01 ) neq = 5.786E-01 ; else if(energy < 7.050E-01 ) neq = 6.019E-01 ; else if(energy < 7.400E-01 ) neq = 6.724E-01 ; else if(energy < 7.800E-01 ) neq = 9.209E-01 ; else if(energy < 8.200E-01 ) neq = 1.459E+00 ; else if(energy < 8.600E-01 ) neq = 8.318E-01 ; else if(energy < 9.000E-01 ) neq = 9.434E-01 ; else if(energy < 9.400E-01 ) neq = 1.172E+00 ; else if(energy < 9.800E-01 ) neq = 1.174E+00 ; else if(energy < 1.050E+00 ) neq = 8.020E-01 ; else if(energy < 1.150E+00 ) neq = 6.578E-01 ; else if(energy < 1.250E+00 ) neq = 9.680E-01 ; else if(energy < 1.350E+00 ) neq = 9.410E-01 ; else if(energy < 1.450E+00 ) neq = 1.079E+00 ; else if(energy < 1.550E+00 ) neq = 1.128E+00 ; else if(energy < 1.650E+00 ) neq = 1.766E+00 ; else if(energy < 1.750E+00 ) neq = 8.366E-01 ; else if(energy < 1.850E+00 ) neq = 1.411E+00 ; else if(energy < 1.950E+00 ) neq = 1.393E+00 ; else if(energy < 2.050E+00 ) neq = 1.022E+00 ; else if(energy < 2.150E+00 ) neq = 1.159E+00 ; else if(energy < 2.250E+00 ) neq = 1.126E+00 ; else if(energy < 2.350E+00 ) neq = 1.106E+00 ; else if(energy < 2.450E+00 ) neq = 1.267E+00 ; else if(energy < 2.550E+00 ) neq = 1.384E+00 ; else if(energy < 2.650E+00 ) neq = 1.238E+00 ; else if(energy < 2.750E+00 ) neq = 1.153E+00 ; else if(energy < 2.850E+00 ) neq = 1.437E+00 ; else if(energy < 2.950E+00 ) neq = 1.061E+00 ; else if(energy < 3.050E+00 ) neq = 1.278E+00 ; else if(energy < 3.150E+00 ) neq = 1.427E+00 ; else if(energy < 3.250E+00 ) neq = 1.281E+00 ; else if(energy < 3.350E+00 ) neq = 1.217E+00 ; else if(energy < 3.450E+00 ) neq = 1.259E+00 ; else if(energy < 3.550E+00 ) neq = 1.213E+00 ; else if(energy < 3.650E+00 ) neq = 7.506E-01 ; else if(energy < 3.750E+00 ) neq = 1.235E+00 ; else if(energy < 3.850E+00 ) neq = 1.187E+00 ; else if(energy < 3.950E+00 ) neq = 1.452E+00 ; else if(energy < 4.050E+00 ) neq = 1.431E+00 ; else if(energy < 4.150E+00 ) neq = 1.147E+00 ; else if(energy < 4.250E+00 ) neq = 1.782E+00 ; else if(energy < 4.350E+00 ) neq = 1.442E+00 ; else if(energy < 4.450E+00 ) neq = 1.509E+00 ; else if(energy < 4.550E+00 ) neq = 1.498E+00 ; else if(energy < 4.650E+00 ) neq = 1.673E+00 ; else if(energy < 4.750E+00 ) neq = 2.003E+00 ; else if(energy < 4.850E+00 ) neq = 1.717E+00 ; else if(energy < 4.950E+00 ) neq = 1.579E+00 ; else if(energy < 5.050E+00 ) neq = 1.612E+00 ; else if(energy < 5.150E+00 ) neq = 1.831E+00 ; else if(energy < 5.250E+00 ) neq = 1.605E+00 ; else if(energy < 5.350E+00 ) neq = 1.336E+00 ; else if(energy < 5.450E+00 ) neq = 1.277E+00 ; else if(energy < 5.550E+00 ) neq = 1.567E+00 ; else if(energy < 5.650E+00 ) neq = 1.616E+00 ; else if(energy < 5.750E+00 ) neq = 1.929E+00 ; else if(energy < 5.850E+00 ) neq = 1.794E+00 ; else if(energy < 5.950E+00 ) neq = 1.468E+00 ; else if(energy < 6.050E+00 ) neq = 1.666E+00 ; else if(energy < 6.150E+00 ) neq = 1.364E+00 ; else if(energy < 6.250E+00 ) neq = 1.900E+00 ; else if(energy < 6.350E+00 ) neq = 1.641E+00 ; else if(energy < 6.450E+00 ) neq = 1.514E+00 ; else if(energy < 6.550E+00 ) neq = 1.314E+00 ; else if(energy < 6.650E+00 ) neq = 1.611E+00 ; else if(energy < 6.750E+00 ) neq = 1.778E+00 ; else if(energy < 6.850E+00 ) neq = 1.577E+00 ; else if(energy < 6.950E+00 ) neq = 1.521E+00 ; else if(energy < 7.050E+00 ) neq = 1.770E+00 ; else if(energy < 7.150E+00 ) neq = 1.470E+00 ; else if(energy < 7.250E+00 ) neq = 1.803E+00 ; else if(energy < 7.350E+00 ) neq = 1.818E+00 ; else if(energy < 7.450E+00 ) neq = 1.793E+00 ; else if(energy < 7.550E+00 ) neq = 1.757E+00 ; else if(energy < 7.650E+00 ) neq = 1.794E+00 ; else if(energy < 7.750E+00 ) neq = 1.866E+00 ; else if(energy < 7.850E+00 ) neq = 1.861E+00 ; else if(energy < 7.950E+00 ) neq = 1.992E+00 ; else if(energy < 8.050E+00 ) neq = 1.806E+00 ; else if(energy < 8.150E+00 ) neq = 1.663E+00 ; else if(energy < 8.250E+00 ) neq = 1.743E+00 ; else if(energy < 8.350E+00 ) neq = 1.768E+00 ; else if(energy < 8.450E+00 ) neq = 1.782E+00 ; else if(energy < 8.550E+00 ) neq = 1.771E+00 ; else if(energy < 8.650E+00 ) neq = 1.769E+00 ; else if(energy < 8.750E+00 ) neq = 1.589E+00 ; else if(energy < 8.850E+00 ) neq = 1.685E+00 ; else if(energy < 8.950E+00 ) neq = 1.871E+00 ; else if(energy < 9.050E+00 ) neq = 1.876E+00 ; else if(energy < 9.150E+00 ) neq = 1.686E+00 ; else if(energy < 9.250E+00 ) neq = 1.619E+00 ; else if(energy < 9.350E+00 ) neq = 1.751E+00 ; else if(energy < 9.450E+00 ) neq = 1.824E+00 ; else if(energy < 9.550E+00 ) neq = 1.767E+00 ; else if(energy < 9.650E+00 ) neq = 1.669E+00 ; else if(energy < 9.750E+00 ) neq = 1.707E+00 ; else if(energy < 9.850E+00 ) neq = 1.754E+00 ; else if(energy < 9.950E+00 ) neq = 1.780E+00 ; else if(energy < 1.005E+01 ) neq = 1.760E+00 ; else if(energy < 1.015E+01 ) neq = 1.782E+00 ; else if(energy < 1.025E+01 ) neq = 1.764E+00 ; else if(energy < 1.035E+01 ) neq = 1.772E+00 ; else if(energy < 1.045E+01 ) neq = 1.709E+00 ; else if(energy < 1.055E+01 ) neq = 1.696E+00 ; else if(energy < 1.065E+01 ) neq = 1.674E+00 ; else if(energy < 1.075E+01 ) neq = 1.726E+00 ; else if(energy < 1.085E+01 ) neq = 1.765E+00 ; else if(energy < 1.095E+01 ) neq = 1.765E+00 ; else if(energy < 1.105E+01 ) neq = 1.736E+00 ; else if(energy < 1.115E+01 ) neq = 1.680E+00 ; else if(energy < 1.125E+01 ) neq = 1.692E+00 ; else if(energy < 1.135E+01 ) neq = 1.735E+00 ; else if(energy < 1.145E+01 ) neq = 1.765E+00 ; else if(energy < 1.155E+01 ) neq = 1.778E+00 ; else if(energy < 1.165E+01 ) neq = 1.767E+00 ; else if(energy < 1.175E+01 ) neq = 1.751E+00 ; else if(energy < 1.185E+01 ) neq = 1.754E+00 ; else if(energy < 1.195E+01 ) neq = 1.767E+00 ; else if(energy < 1.205E+01 ) neq = 1.784E+00 ; else if(energy < 1.215E+01 ) neq = 1.812E+00 ; else if(energy < 1.225E+01 ) neq = 1.816E+00 ; else if(energy < 1.235E+01 ) neq = 1.822E+00 ; else if(energy < 1.245E+01 ) neq = 1.791E+00 ; else if(energy < 1.255E+01 ) neq = 1.757E+00 ; else if(energy < 1.265E+01 ) neq = 1.837E+00 ; else if(energy < 1.275E+01 ) neq = 1.856E+00 ; else if(energy < 1.285E+01 ) neq = 1.850E+00 ; else if(energy < 1.295E+01 ) neq = 1.853E+00 ; else if(energy < 1.305E+01 ) neq = 1.862E+00 ; else if(energy < 1.315E+01 ) neq = 1.828E+00 ; else if(energy < 1.325E+01 ) neq = 1.793E+00 ; else if(energy < 1.335E+01 ) neq = 1.806E+00 ; else if(energy < 1.345E+01 ) neq = 1.839E+00 ; else if(energy < 1.355E+01 ) neq = 1.802E+00 ; else if(energy < 1.365E+01 ) neq = 1.780E+00 ; else if(energy < 1.375E+01 ) neq = 1.804E+00 ; else if(energy < 1.385E+01 ) neq = 1.870E+00 ; else if(energy < 1.395E+01 ) neq = 1.808E+00 ; else if(energy < 1.405E+01 ) neq = 1.787E+00 ; else if(energy < 1.415E+01 ) neq = 1.817E+00 ; else if(energy < 1.425E+01 ) neq = 1.830E+00 ; else if(energy < 1.435E+01 ) neq = 1.843E+00 ; else if(energy < 1.445E+01 ) neq = 1.845E+00 ; else if(energy < 1.455E+01 ) neq = 1.791E+00 ; else if(energy < 1.465E+01 ) neq = 1.765E+00 ; else if(energy < 1.475E+01 ) neq = 1.774E+00 ; else if(energy < 1.485E+01 ) neq = 1.808E+00 ; else if(energy < 1.495E+01 ) neq = 1.813E+00 ; else if(energy < 1.505E+01 ) neq = 1.773E+00 ; else if(energy < 1.515E+01 ) neq = 1.801E+00 ; else if(energy < 1.525E+01 ) neq = 1.842E+00 ; else if(energy < 1.535E+01 ) neq = 1.841E+00 ; else if(energy < 1.545E+01 ) neq = 1.828E+00 ; else if(energy < 1.555E+01 ) neq = 1.808E+00 ; else if(energy < 1.565E+01 ) neq = 1.806E+00 ; else if(energy < 1.575E+01 ) neq = 1.817E+00 ; else if(energy < 1.585E+01 ) neq = 1.858E+00 ; else if(energy < 1.595E+01 ) neq = 1.891E+00 ; else if(energy < 1.605E+01 ) neq = 1.899E+00 ; else if(energy < 1.615E+01 ) neq = 1.899E+00 ; else if(energy < 1.625E+01 ) neq = 1.893E+00 ; else if(energy < 1.635E+01 ) neq = 1.880E+00 ; else if(energy < 1.645E+01 ) neq = 1.869E+00 ; else if(energy < 1.655E+01 ) neq = 1.918E+00 ; else if(energy < 1.665E+01 ) neq = 1.970E+00 ; else if(energy < 1.675E+01 ) neq = 1.952E+00 ; else if(energy < 1.685E+01 ) neq = 1.929E+00 ; else if(energy < 1.695E+01 ) neq = 1.935E+00 ; else if(energy < 1.705E+01 ) neq = 1.940E+00 ; else if(energy < 1.715E+01 ) neq = 1.937E+00 ; else if(energy < 1.725E+01 ) neq = 1.934E+00 ; else if(energy < 1.735E+01 ) neq = 1.935E+00 ; else if(energy < 1.745E+01 ) neq = 1.932E+00 ; else if(energy < 1.755E+01 ) neq = 1.917E+00 ; else if(energy < 1.765E+01 ) neq = 1.906E+00 ; else if(energy < 1.775E+01 ) neq = 1.914E+00 ; else if(energy < 1.785E+01 ) neq = 1.926E+00 ; else if(energy < 1.795E+01 ) neq = 1.951E+00 ; else if(energy < 1.805E+01 ) neq = 1.974E+00 ; else if(energy < 1.815E+01 ) neq = 1.965E+00 ; else if(energy < 1.825E+01 ) neq = 1.950E+00 ; else if(energy < 1.835E+01 ) neq = 1.938E+00 ; else if(energy < 1.845E+01 ) neq = 1.926E+00 ; else if(energy < 1.855E+01 ) neq = 1.945E+00 ; else if(energy < 1.865E+01 ) neq = 1.986E+00 ; else if(energy < 1.875E+01 ) neq = 2.008E+00 ; else if(energy < 1.885E+01 ) neq = 1.999E+00 ; else if(energy < 1.895E+01 ) neq = 1.989E+00 ; else if(energy < 1.905E+01 ) neq = 1.974E+00 ; else if(energy < 1.915E+01 ) neq = 1.960E+00 ; else if(energy < 1.925E+01 ) neq = 1.955E+00 ; else if(energy < 1.935E+01 ) neq = 1.953E+00 ; else if(energy < 1.945E+01 ) neq = 1.953E+00 ; else if(energy < 1.955E+01 ) neq = 1.955E+00 ; else if(energy < 1.965E+01 ) neq = 1.958E+00 ; else if(energy < 1.975E+01 ) neq = 1.963E+00 ; else if(energy < 1.985E+01 ) neq = 1.969E+00 ; else if(energy < 1.995E+01 ) neq = 1.975E+00 ; return neq; } G4double SteppingAction::GetElectron1MeVNeutronEquivalence(G4double energy) { double neq(0); if(energy< 3.00E-01 ) neq = 3.18E-03 ; else if(energy< 5.00E-01 ) neq = 8.00E-03 ; else if(energy< 7.00E-01 ) neq = 1.14E-02 ; else if(energy< 1.00E+00 ) neq = 1.54E-02 ; else if(energy< 2.00E+00 ) neq = 2.49E-02 ; else if(energy< 3.00E+00 ) neq = 3.13E-02 ; else if(energy< 5.00E+00 ) neq = 3.98E-02 ; else if(energy< 7.00E+00 ) neq = 4.56E-02 ; else if(energy< 1.00E+01 ) neq = 5.16E-02 ; else if(energy< 2.00E+01 ) neq = 6.24E-02 ; else if(energy< 3.00E+01 ) neq = 6.78E-02 ; else if(energy< 5.00E+01 ) neq = 7.34E-02 ; else if(energy< 7.00E+01 ) neq = 7.63E-02 ; else if(energy< 1.00E+02 ) neq = 7.87E-02 ; else if(energy< 2.00E+02 ) neq = 8.11E-02 ; return neq; }