From 6ac128f92793a858444f484cd20561dc1e25861d Mon Sep 17 00:00:00 2001
From: itayhorin <itayhorin@gmail.com>
Date: Tue, 25 Oct 2022 13:42:25 +0300
Subject: [PATCH] implementd two solution for comCM

---
 bla                                    |   0
 include/openmc/tallies/tally_scoring.h |   2 +-
 src/tallies/tally_scoring.cpp          | 283 +++++++++++++++++++++++--
 3 files changed, 272 insertions(+), 13 deletions(-)
 create mode 100644 bla

diff --git a/bla b/bla
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index 1765c333850..888539089f8 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -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.
 //
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 1c3fa852d69..0e727fec1aa 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -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());
@@ -2449,14 +2447,17 @@ 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;
@@ -2464,9 +2465,9 @@ void score_point_tally(Particle& p)
   //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);
@@ -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);
@@ -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 )
@@ -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