Skip to content

Commit

Permalink
implementd two solution for comCM
Browse files Browse the repository at this point in the history
  • Loading branch information
@Itay-max
Itay-max committed Oct 25, 2022
1 parent 23a7d14 commit 6ac128f
Show file tree
Hide file tree
Showing 3 changed files with 272 additions and 13 deletions.
Empty file added bla
View file
Empty file.
2 changes: 1 addition & 1 deletion include/openmc/tallies/tally_scoring.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ void Vcros(double A[4],double B[4],double C[4]);
void Vunit(double A[4] ,double B[4]);


void getMu_COM(double x_det , double y_det , double z_det ,Particle p_col , double awr , double incoming_mass, double ReturnArray[],int diff_mode,double dl );
void getMu_COM(double x_det , double y_det , double z_det ,Particle p_col , double awr , double ReturnArray[],int diff_mode,double dl );

//! Analog tallies are triggered at every collision, not every event.
//
Expand Down
283 changes: 271 additions & 12 deletions src/tallies/tally_scoring.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2432,14 +2432,12 @@ void score_point_tally(Particle& p)
std::cout << "mass in ev " << p.getMass() << std::endl ;
// Determine the collision estimate of the flux
bool verbose=true;
double ReturnArray[2];
double ReturnArray[4];
double flux = 0.0;
const auto& nuc {data::nuclides[p.event_nuclide()]};
double awr = nuc->awr_;
//auto [value1, value2] = getMu_COM(0,0,0,p,awr,p.getMass(),ReturnArray);
//double ResultsArray = getMu_COM(0,0,0,p,awr,p.getMass());
double dl = 1e-5;
getMu_COM(0,0,0,p,awr,p.getMass(),ReturnArray , 0, dl);
getMu_COM(0,0,0,p,awr,ReturnArray , 0, dl);
const auto& rx {nuc->reactions_[0]};
auto& d = rx->products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
Expand All @@ -2449,24 +2447,27 @@ void score_point_tally(Particle& p)
double theta_lab = std::acos(cos_lab);
double sin_lab = std::sin(theta_lab);
double mu_COM_or_itay = (std::sqrt(awr*awr - sin_lab*sin_lab)*cos_lab - sin_lab*sin_lab)/awr;
double mu_COM = ReturnArray[0];
double ReturnArrayPlus[2];
double ReturnArrayMinus[2];
double mu_COM1 = ReturnArray[0];
double E1 = ReturnArray[2];
double mu_COM2 = ReturnArray[1];
double E2 = ReturnArray[3];
//double ReturnArrayPlus[4];
//double ReturnArrayMinus[4];
//getMu_COM(0,0,0,p,awr,p.getMass(),ReturnArrayPlus ,1 , 1e-5 );
//getMu_COM(0,0,0,p,awr,p.getMass(),ReturnArrayMinus ,-1 , 1e-5 );
//double MuPlus = ReturnArrayPlus[0]; // not sure about changing mu lab correctly
//double MuMinus = ReturnArrayMinus[0];
double theta_pdf = d_->angle().get_pdf_value(p.E_last(),mu_COM,p.current_seed());
double theta_pdf = d_->angle().get_pdf_value(p.E_last(),mu_COM1,p.current_seed());

//double derivative = std::abs(MuPlus-MuMinus)/(2*dl);
double derivative =1;
double theta_pdf_lab = theta_pdf * derivative;
//double E_ghost = p.E_last()*(1+awr*awr+2*awr*mu_COM)/(1+awr)/(1+awr);

//double m_incoming =MASS_NEUTRON_EV;
double E_ghost = ReturnArray[1];
double E_ghost = E1;

if(!std::isnan(mu_COM))
if(true)//(!std::isnan(mu_COM1))
{
Particle ghost_particle=Particle();
ghost_particle.initilze_ghost_particle(p,u_lab,E_ghost);
Expand Down Expand Up @@ -2519,7 +2520,8 @@ void score_point_tally(Particle& p)
fmt::print("parent particle E = {}\n",p.E_last());
fmt::print("ghost particle E = {}\n",ghost_particle.E());
fmt::print("ghost particle u = {0} , {1} , {2}\n",ghost_particle.u().x,ghost_particle.u().y,ghost_particle.u().z);
fmt::print("ghost particle Mu COM Arik= {}\n",mu_COM);
fmt::print("ghost particle Mu COM 1 Arik= {}\n",mu_COM1);
fmt::print("ghost particle Mu COM 2 Arik= {}\n",mu_COM2);
fmt::print("ghost particle Mu COM Or and Itay = {}\n",mu_COM_or_itay);
fmt::print("flux = {}\n",flux);
fmt::print("theta cm pdf ={} \n",theta_pdf);
Expand Down Expand Up @@ -2635,7 +2637,7 @@ void score_surface_tally(Particle& p, const vector<int>& tallies)
}



/*
void getMu_COM(double x_det , double y_det , double z_det ,Particle p_col , double awr , double incoming_mass ,double ReturnArray[],int diff_mode,double dl )
Expand Down Expand Up @@ -2882,6 +2884,263 @@ void Vcros(double A[4],double B[4],double C[4])
return;
}
*/


//Author: Arik Kreisel


void getMu_COM(double x_det , double y_det , double z_det ,Particle p_col , double awr ,double ReturnArray[],int diff_mode,double dl )
{

double cs=0;
int m=0;

double m1= p_col.getMass()/1e6; // mass of incoming particle in MeV
double E1_tot = p_col.E_last()/1e6 + m1; // total Energy of incoming particle in MeV
double p1_tot = std::sqrt(E1_tot*E1_tot - m1*m1);
double r1[4] = {0, x_det, y_det, z_det}; // detector position lab {ignor, x, y, z}
double r2[4]= {0, p_col.r().x, p_col.r().y, p_col.r().z}; // collision position lab {ignor, x, y, z}
double r3[4]; // r1-r2 vector from collision to detector
double m2= m1*awr; // mass of target matirial
double m3= m1; // mass of outgoing particle to detector
double m4= m2; // mass of recoil target system
double p1[3]={p1_tot*p_col.u_last().x,p1_tot* p_col.u_last().y,p1_tot* p_col.u_last().z}; // 3 momentum of incoming particle
double p2[3]={0, 0, 0}; //3 momentum of target in lab

// calculate
double Fp1[4]; //four momentum of incoming particle in LAB
double Fp2[4]; //four momentum of target in LAB
double Fp3[4]; //four momentum of particle going to the detector in LAB
double UFp3[4]; //unit 3 momentum of particle going to the detector
double CM[4]; //four momentum of center of mass frame in lab
double LAB[4]; //four momentum of lab in center of mass frame
double mCM; // mass of center of mass system
double pCM; // momentum of center of mass system
double p3LAB[2]; // momentum of out particle

double CMFp1[4]; //four momentum of incoming particle in CM
double CMFp2[4]; //four momentum of incoming target in CM
double CMFp3[4]; //four momentum of out particle in CM
double CMFp4[4]; //four momentum of out target in CM
double Fp4[4]; //four momentum of out target in LAB
double CMUp1[4]; //unit three vector of incoming particle in CM
double CMUp3[4]; //unit three vector of out particle in CM
double cosCM; // cosine of out going particle to detector in CM frame
double cosLAB; // cosine of out going particle to detector in LAB frame
double CME3; // Energy of out particle in CM
double CMp3; // momentum of out particle in CM
double Ur3[4], UCM[4], ULAB[4];
double aa, bb, cc;


Fp1[0]=0;
Fp2[0]=0;
CM[0]=0;
LAB[0]=0;
r3[0]=0;


for(int i=0; i<3; i++){
Fp1[i+1]=p1[i];
Fp2[i+1]=p2[i];
CM[i+1]=Fp1[i+1]+Fp2[i+1];
LAB[i+1]=-CM[i+1];
r3[i+1]=r1[i+1]-r2[i+1];
}

Fp1[0]=sqrt(Fp1[1]*Fp1[1]+Fp1[2]*Fp1[2]+Fp1[3]*Fp1[3]+m1*m1);
Fp2[0]=sqrt(Fp2[1]*Fp2[1]+Fp2[2]*Fp2[2]+Fp2[3]*Fp2[3]+m2*m2);
CM[0]=Fp1[0]+Fp2[0];
LAB[0]=CM[0];
r3[0]=0;

mCM=sqrt(CM[0]*CM[0]-CM[1]*CM[1]-CM[2]*CM[2]-CM[3]*CM[3]);
pCM=sqrt(CM[1]*CM[1]+CM[2]*CM[2]+CM[3]*CM[3]);
CME3=(mCM*mCM-m4*m4+m3*m3)/(2*mCM); // energy of out going particle in CM
CMp3=sqrt(CME3*CME3-m3*m3);

std::cout<<" mCM= "<<mCM<<" pCM= "<<pCM<<" CME3= "<<CME3<<std::endl;

boostf(CM,Fp1,CMFp1);
boostf(CM,Fp2,CMFp2);

Vunit(CM,UCM);
Vunit(LAB,ULAB);
Vunit(r3,Ur3);
cosLAB=Vdot(UCM,Ur3);

std::cout<<"cosLAB= "<<cosLAB<<std::endl;
Vunit(CMFp1,CMUp1);


aa=pCM*pCM*cosLAB*cosLAB-CM[0]*CM[0];
bb=2*pCM*cosLAB*CME3*mCM;
cc=CME3*mCM*CME3*mCM-m3*m3*CM[0]*CM[0];

int j=1;

if(bb*bb-4*aa*cc<0) {
std::cout<<" detector out of range" <<std::endl;
return; //continue;
}
p3LAB[0]=(-bb+sqrt(bb*bb-4*aa*cc))/2.0/aa;
if(p3LAB[0]<=0) {p3LAB[0]=(-bb-sqrt(bb*bb-4*aa*cc))/2/aa;
if(p3LAB[0]<=0) {
std::cout<<" detector out of range" <<std::endl;
return; //continue;
}
}
else if(p3LAB[0]>0){p3LAB[1]=(-bb-sqrt(bb*bb-4*aa*cc))/2/aa;
if(p3LAB[1]>0) j=j+1;
}

std::cout<<" p3LAB1= "<<(-bb+sqrt(bb*bb-4*aa*cc))/2.0/aa<<" p3LAB2= "<<(-bb-sqrt(bb*bb-4*aa*cc))/2.0/aa<<std::endl;

for (int l=0;l<j;l++){

std::cout<<"l= "<<l<<std::endl;

Fp3[0]=sqrt(p3LAB[l]*p3LAB[l]+m3*m3);
for(int i=0; i<3; i++){
Fp3[i+1]=p3LAB[l]*Ur3[i+1];
}

boostf(CM,Fp3,CMFp3);
Vunit(CMFp3,CMUp3);
cosCM=Vdot(UCM,CMUp3); // input to openMC Cross section calculation
ReturnArray[l] = cosCM;
ReturnArray[l+2] = (Fp3[0]-m3)*1e6; //retrun Energy in eV

std::cout<<"cosCM= "<<cosCM<<std::endl;

std::cout<<" ----------- sanity tests for solution "<< l <<" -----------------------------------"<<std::endl;

Vunit(Fp3,UFp3);

std::cout<<Vdot(UCM,UFp3)<<" = cosLAB = "<<cosLAB<<std::endl;

std::cout<<CMFp3[0]<<" m3 energy in CM "<<CME3<<std::endl;

CMFp4[0]=0;
for(int i=0; i<3; i++){
CMFp4[i+1]=-CMFp3[i+1];
}
CMFp4[0]=sqrt(m4*m4+CMFp4[1]*CMFp4[1]+CMFp4[2]*CMFp4[2]+CMFp4[3]*CMFp4[3]);
boostf(LAB,CMFp4,Fp4);

for(int i=0; i<4; i++){
std::cout<<i<<" "<<CM[i]<<" = energy momentum conserv in LAB = "<<Fp4[i]+Fp3[i]<<std::endl;
std::cout<<i<<" "<<CMFp1[i]+CMFp2[i]<<" = energy momentum conserv in CM = "<<CMFp4[i]+CMFp3[i]<<std::endl;
}



double test[4];
boostf(LAB,CMFp3,test);
for (int i=0; i<4; i++){
std::cout<<i<<" "<<test[i]<<std::endl;
std::cout<<i<<" "<<Fp3[i]<<std::endl;
}

std::cout <<" m3 enerrgy in lab ="<<Fp3[0]<<std::endl;
std::cout <<" m4 enerrgy in lab ="<<Fp4[0]<<std::endl;
std::cout <<" sqrt(E3^2-P3^2) in lab ="<<sqrt(Fp3[0]*Fp3[0]-Fp3[1]*Fp3[1]-Fp3[2]*Fp3[2]-Fp3[3]*Fp3[3])<<std::endl;
std::cout <<" sqrt(E4^2-P4^2) in lab ="<<sqrt(Fp4[0]*Fp4[0]-Fp4[1]*Fp4[1]-Fp4[2]*Fp4[2]-Fp4[3]*Fp4[3])<<std::endl;
std::cout <<" m1 enerrgy in lab ="<<Fp1[0]<<std::endl;

double tanLab=sqrt(1-cosCM*cosCM)/(CM[0]/mCM*(pCM*CME3/CM[0]/CMp3+cosCM));
double fCOSlab=cos(atan(tanLab));


std::cout<<"cos( atan( tanLab= "<<cos(atan(tanLab))<<std::endl;

}

}


void boostf( double A[4], double B[4], double X[4])
{
//
// boosts B(labfram) to A rest frame and gives output in X
//
double W;
int j;

if ((A[0]*A[0]-A[1]*A[1]-A[2]*A[2]-A[3]*A[3])<=0) {
std::cout <<"negative sqrt in boostf"<<A[0]<<A[1]<<A[2]<<A[3]<<std::endl;}

W=sqrt(A[0]*A[0]-A[1]*A[1]-A[2]*A[2]-A[3]*A[3]);

if(W==0 || W==(-A[0])) std::cout <<"divid by 0 in boostf"<<std::endl;

X[0]=(B[0]*A[0]-B[3]*A[3]-B[2]*A[2]-B[1]*A[1])/W;
for(j=1; j<=3; j++) {
X[j]=B[j]-A[j]*(B[0]+X[0])/(A[0]+W);
}

return;
}




double Vdot(double A[4],double B[4])
{
int j;

double dot = 0;

for(j=1; j<=3; j++) {
dot = dot + A[j]*B[j];
}


return dot;
}


void Vunit(double A[4] ,double B[4])
{
double fff;
int j;

fff = 0;

for(j=1; j<=3; j++) {
fff = fff + A[j]*A[j];
}

if (fff==0) {
std::cout <<"in vunit divid by zero" << std::endl;
return;
}

for(j=1; j<=3; j++) {
B[j] = A[j]/sqrt(fff);
}
B[0] = 0;

return;
}


void Vcros(double A[4],double B[4],double C[4])
{
C[1] = A[2]*B[3]-A[3]*B[2];
C[2] = A[3]*B[1]-A[1]*B[3];
C[3] = A[1]*B[2]-A[2]*B[1];
C[0] = 0;

if (C[1]==0 && C[2]==0 && C[3]==0.) {
std::cout << "vcross zero" << std::endl;
}

return;
}




} // namespace openmc

0 comments on commit 6ac128f

Please sign in to comment.