Skip to content

geant
geant
Commit 398cd2ce authored by Mihaly Novak's avatar Mihaly Novak
Browse files
Options

Merge branch

parents 13c337a0 cc90bae4
Hide whitespace changes
Inline Side-by-side
Showing with 115 additions and 113 deletions
test-small/src/TabulatedDataManager.cc
View file @ 398cd2ce
......@@ -452,9 +452,9 @@ Int_t TabulatedDataManager::SampleInteraction( const G4int imat,
// kinEnergy = atrack.GetKineticEnergy()/CLHEP::GeV; // from MeV->GeV
// TO BE THE SAME AS G4 WITHOUT integral approach !!! IT IS DIFFERENT IT PROTO.
// BUT NO EFFECT IF THE TRACKING LIMIT IS HIGHER THAN 3keV
if(fgIsUseRange)
kinEnergy = atrack.GetStep()->GetPreStepPoint()->GetKineticEnergy()/GeV;
else
// if(fgIsUseRange)
// kinEnergy = atrack.GetStep()->GetPreStepPoint()->GetKineticEnergy()/GeV;
// else
kinEnergy = atrack.GetStep()->GetPostStepPoint()->GetKineticEnergy()/GeV;
// sampling element for intercation based on element-wise relative tot-xsecs
......
xsec/src/TMXsec.cxx
View file @ 398cd2ce
......@@ -277,11 +277,11 @@ Float_t TMXsec::Xlength(Int_t part, Float_t en, Double_t ptot) {
// Double_t en2 = en1*fEDelta;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("Xlength","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return TMath::Limits<Float_t>::Max();
}
// if(en1>en || en2<en) {
// Error("Xlength","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return TMath::Limits<Float_t>::Max();
// }
Double_t xrat = (en2-en)/(en2-en1);
Double_t x = xrat*fTotXL[part*fNEbins+ibin]+(1-xrat)*fTotXL[part*fNEbins+ibin+1];
return x;
......@@ -419,11 +419,11 @@ Float_t TMXsec::DEdx(Int_t part, Float_t en) {
// Double_t en2 = en1*fEDelta;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("DEdx","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return TMath::Limits<Float_t>::Max();
}
// if(en1>en || en2<en) {
// Error("DEdx","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return TMath::Limits<Float_t>::Max();
// }
Double_t xrat = (en2-en)/(en2-en1);
return xrat*fDEdx[part*fNEbins+ibin]+(1-xrat)*fDEdx[part*fNEbins+ibin+1];
......@@ -475,11 +475,11 @@ Float_t TMXsec::Range(Int_t part, Float_t en) {
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("Range","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return -1.0;
}
// if(en1>en || en2<en) {
// Error("Range","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return -1.0;
// }
Double_t xrat = (en2-en)/(en2-en1);
return xrat*fRange[part*fNEbins+ibin]+(1-xrat)*fRange[part*fNEbins+ibin+1];
......@@ -619,12 +619,12 @@ TEXsec* TMXsec::SampleInt(Int_t part, Double_t en, Int_t &reac, Double_t ptot) {
// Double_t en2 = en1*fEDelta;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SampleInt","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
reac=-1;
return 0;
}
// if(en1>en || en2<en) {
// Error("SampleInt","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// reac=-1;
// return 0;
// }
Int_t iel=-1;
if(fNElems==1) {
iel=0;
......
xsec/src/TPFstate.cxx
View file @ 398cd2ce
......@@ -153,11 +153,11 @@ Bool_t TPFstate::SampleReac(Int_t preac, Float_t en, Int_t& npart, Float_t& weig
ibin = ibin<fNEbins-1?ibin:fNEbins-2;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return kFALSE;
}
// if(en1>en || en2<en) {
// Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return kFALSE;
// }
Double_t xrat = (en2-en)/(en2-en1);
if(eta>xrat) ++ibin;
ebinindx = ibin;
......@@ -191,11 +191,11 @@ Bool_t TPFstate::SampleReac(Int_t preac, Float_t en, Int_t& npart, Float_t& weig
ibin = ibin<fNEbins-1?ibin:fNEbins-2;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return kFALSE;
}
// if(en1>en || en2<en) {
// Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return kFALSE;
// }
Double_t xrat = (en2-en)/(en2-en1);
if(randn1>xrat) ++ibin;
Int_t ipoint = rnumber*fNEbins + ibin;
......@@ -258,11 +258,11 @@ Bool_t TPFstate::GetReac(Int_t preac, Float_t en, Int_t ifs, Int_t& npart, Float
ibin = ibin<fNEbins-1?ibin:fNEbins-2;
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return kFALSE;
}
// if(en1>en || en2<en) {
// Error("SetFinState","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return kFALSE;
// }
if(en-en1>en2-en) ++ibin;
Int_t ipoint = rnumber*fNEbins + ibin;
// in case of any problems with the fstate sampling the primary will be
......@@ -332,10 +332,10 @@ Bool_t TPFstate::Resample() {
ibin = ibin<oNEbins-1?ibin:oNEbins-2;
Double_t en1 = oGrid[ibin];
Double_t en2 = oGrid[ibin+1];
if(en1>en || en<en) {
Error("Interp","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
}
// if(en1>en || en<en) {
// Error("Interp","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// }
Double_t xrat = (en2-en)/(en2-en1);
if(xrat<0.5) obins[ien]=ibin;
else obins[ien]=ibin+1;
......
xsec/src/TPXsec.cxx
View file @ 398cd2ce
......@@ -93,10 +93,10 @@ void TPXsec::Interp(Double_t egrid[], Float_t value[], Int_t nbins,
ibin = ibin<nbins-1?ibin:nbins-2;
Double_t en1 = egrid[ibin];
Double_t en2 = egrid[ibin+1];
if(en1>en || en2<en) {
Error("Interp","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
}
// if(en1>en || en2<en) {
// Error("Interp","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// }
Double_t xrat = (en2-en)/(en2-en1);
for(Int_t ival=0; ival<stride; ++ival) {
result[ival] = xrat*value[ibin*stride+ival]+
......@@ -289,11 +289,11 @@ Float_t TPXsec::DEdx(Double_t en) const {
// Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return 0;
}
// if(en1>en || en2<en) {
// Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return 0;
// }
Double_t xrat = (en2-en)/(en2-en1);
Double_t dedx = xrat*fdEdx[ibin]+(1-xrat)*fdEdx[ibin+1];
/* printf("ibin %d en1 %f en %f en2 %f xs1 %f xs2 %f xrat %f xsec %f\n",
......@@ -316,11 +316,11 @@ Bool_t TPXsec::MS(Double_t en, Float_t &ang, Float_t &asig,
// Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return 0;
}
// if(en1>en || en2<en) {
// Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return 0;
// }
Double_t xrat = (en2-en)/(en2-en1);
ang = xrat*fMSangle[ibin]+(1-xrat)*fMSangle[ibin+1];
asig = xrat*fMSansig[ibin]+(1-xrat)*fMSansig[ibin+1];
......@@ -339,11 +339,12 @@ Int_t TPXsec::SampleReac(Double_t en) const {
//Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SampleReac","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return 0;
}
// if(en1>en || en2<en) {
// Error("SampleReac","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return -1;
// }
Double_t xrat = (en2-en)/(en2-en1);
Double_t xnorm = 1.;
while(1) {
......@@ -367,11 +368,12 @@ Int_t TPXsec::SampleReac(Double_t en, Double_t randn) const {
//Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("SampleReac","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return 0;
}
// if(en1>en || en2<en) {
// Error("SampleReac","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return -1;
// }
Double_t xrat = (en2-en)/(en2-en1);
Double_t xnorm = 1.;
while(1) {
......@@ -399,11 +401,11 @@ Bool_t TPXsec::XS_v(Int_t npart, Int_t rindex, const Double_t en[], Double_t lam
// Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>ene || en2<ene) {
Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, ene, en2);
return 0;
}
// if(en1>ene || en2<ene) {
// Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, ene, en2);
// return 0;
// }
Double_t xrat = (en2-ene)/(en2-en1);
Double_t xtot = (xrat*fTotXs[ibin]+(1-xrat)*fTotXs[ibin+1]);
......@@ -434,11 +436,11 @@ Float_t TPXsec::XS(Int_t rindex, Double_t en) const {
// Double_t en2 = fEmin*TMath::Exp((ibin+1)/fEilDelta);
Double_t en1 = fEGrid[ibin];
Double_t en2 = fEGrid[ibin+1];
if(en1>en || en2<en) {
Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
ibin, en1, en, en2);
return 0;
}
// if(en1>en || en2<en) {
// Error("XS","Wrong bin %d in interpolation: should be %f < %f < %f\n",
// ibin, en1, en, en2);
// return 0;
// }
Double_t xrat = (en2-en)/(en2-en1);
Double_t xtot = (xrat*fTotXs[ibin]+(1-xrat)*fTotXs[ibin+1]);
......
xsec/src/TTabPhysMgr.cxx
View file @ 398cd2ce
......@@ -607,58 +607,58 @@ void TTabPhysMgr::GetRestFinStates(Int_t partindex, TMXsec *mxs,
Double_t randDirX = randSinTheta*TMath::Cos(randPhi);
Double_t randDirY = randSinTheta*TMath::Sin(randPhi);
GeantTrack &gTrack1 = GeantPropagator::Instance()->GetTempTrack(tid);
GeantTrack &gTrack2 = GeantPropagator::Instance()->GetTempTrack(tid);
// need to do it one-by-one
// 1. gamma
GeantTrack &gTrack1 = propagator->GetTempTrack(tid);
//set the new track properties: 2 gamma with m_{e}*c*c
gTrack1.fEvent = gTrack2.fEvent = tracks.fEventV[iintrack];
gTrack1.fEvslot = gTrack2.fEvslot = tracks.fEvslotV[iintrack];
gTrack1.fEvent = tracks.fEventV[iintrack];
gTrack1.fEvslot = tracks.fEvslotV[iintrack];
// gTrack.fParticle = nTotSecPart; //index of this particle
gTrack1.fPDG = gTrack2.fPDG = 22; //gamma PDG code
gTrack1.fG5code = gTrack2.fG5code = TPartIndex::I()->GetSpecGVIndex(2); //gamma G5 index
gTrack1.fEindex = gTrack2.fEindex = 0;
gTrack1.fCharge = gTrack2.fCharge = 0.; // e+ charge
gTrack1.fProcess = gTrack2.fProcess = 0;
gTrack1.fIzero = gTrack2.fIzero = 0;
gTrack1.fNsteps = gTrack2.fNsteps = 0;
gTrack1.fPDG = 22; //gamma PDG code
gTrack1.fG5code = TPartIndex::I()->GetSpecGVIndex(2); //gamma G5 index
gTrack1.fEindex = 0;
gTrack1.fCharge = 0.; // charge
gTrack1.fProcess = 0;
gTrack1.fIzero = 0;
gTrack1.fNsteps = 0;
// gTrack.fSpecies = 0;
gTrack1.fStatus = gTrack2.fStatus = kNew; //status of this particle
gTrack1.fMass = gTrack2.fMass = 0.; //mass of this particle
gTrack1.fXpos = gTrack2.fXpos = tracks.fXposV[iintrack]; //rx of this particle (same as parent)
gTrack1.fYpos = gTrack2.fYpos = tracks.fYposV[iintrack]; //ry of this particle (same as parent)
gTrack1.fZpos = gTrack2.fZpos = tracks.fZposV[iintrack]; //rz of this particle (same as parent)
gTrack1.fStatus = kNew; //status of this particle
gTrack1.fMass = 0.; //mass of this particle
gTrack1.fXpos = tracks.fXposV[iintrack]; //rx of this particle (same as parent)
gTrack1.fYpos = tracks.fYposV[iintrack]; //ry of this particle (same as parent)
gTrack1.fZpos = tracks.fZposV[iintrack]; //rz of this particle (same as parent)
gTrack1.fXdir = randDirX;
gTrack2.fXdir = -1.*randDirX;
gTrack1.fYdir = randDirY;
gTrack2.fYdir = -1.*randDirY;
gTrack1.fZdir = randDirZ;
gTrack2.fZdir = -1.*randDirZ;
gTrack1.fP = gTrack2.fP = mecc; //momentum of this particle
gTrack1.fE = gTrack2.fE = mecc; //total E of this particle
gTrack1.fEdep = gTrack2.fEdep = 0.;
gTrack1.fPstep = gTrack2.fPstep = 0.;
gTrack1.fStep = gTrack2.fStep = 0.;
gTrack1.fSnext = gTrack2.fSnext = 0.;
gTrack1.fSafety = gTrack2.fSafety = tracks.fSafetyV[iintrack];
gTrack1.fFrombdr = gTrack2.fFrombdr= tracks.fFrombdrV[iintrack];
gTrack1.fPending = gTrack2.fPending= kFALSE;
*gTrack1.fPath = *gTrack2.fPath = *tracks.fPathV[iintrack];
*gTrack1.fNextpath = *gTrack2.fNextpath = *tracks.fPathV[iintrack];
gTrack1.fP = mecc; //momentum of this particle
gTrack1.fE = mecc; //total E of this particle
gTrack1.fEdep = 0.;
gTrack1.fPstep = 0.;
gTrack1.fStep = 0.;
gTrack1.fSnext = 0.;
gTrack1.fSafety = tracks.fSafetyV[iintrack];
gTrack1.fFrombdr = tracks.fFrombdrV[iintrack];
gTrack1.fPending = kFALSE;
*gTrack1.fPath = *tracks.fPathV[iintrack];
*gTrack1.fNextpath = *tracks.fPathV[iintrack];
gPropagator->AddTrack(gTrack1);
tracks.AddTrack(gTrack1);
// 2. gamma : everything is the same but the direction
gTrack1.fXdir = -1.*randDirX;
gTrack1.fYdir = -1.*randDirY;
gTrack1.fZdir = -1.*randDirZ;
gPropagator->AddTrack(gTrack1);
tracks.AddTrack(gTrack1);
gPropagator->AddTrack(gTrack2);
tracks.AddTrack(gTrack2);
nTotSecPart+=2;
return;
} else {
return;
}
}
// If the stopped particle doesn't have nuclear capture at-rest then decay it
if(!fHasNCaptureAtRest[partindex]) {
// Decay at-rest
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment