#include "../Rootheaders/Root_macros.h"
#include "TRandom.h"


float radcon=2*3.14159/360;


// REACTION  6Li(n,t)
/*
char *reac="6Li(n,t)";
float m_a=1.0086649; //amu
float m_A=6.0151223; //amu
float m_B=4.0026032; //amu
float m_b= 3.0160493; //amu
//float m_B=3.0160493; //amu
//float m_b= 4.0026032; //amu
 */

// REACTION  10B(n,a)
/*
char *reac="10B(n,a)";
 float m_a=1.0086649; //amu
 float m_A=10.0129370; //amu
 float m_B=7.016004; //amu
 float m_b= 4.002603; //amu
 */
 
float Ex=0.0; //Excitation Energy 478


// REACTION  35Cl(n,p)
char *reac="35Cl(n,p)";
 float m_a=1.0086649; //amu
 float m_A=34.968852; //amu
 float m_B=34.9690322; //amu
 float m_b= 1.0078250; //amu
 


float E_a=3.00; //MeV

float mass_scaler=931.494061; //MeV/c2



void runlist(){
    
    cout<<reac<<endl;
    
    float Q=(m_a+m_A-m_b-m_B)*mass_scaler-Ex;
    cout<< " Q value of reaction ="<<Q<<" MeV"<<endl;
    int i=0;
    
    while(i<=180){
        
        float angle=float(i);
        float theta_cm=angle*2*3.14159/360;
        
        i=i+5;
        
        float term_1=m_a*m_b*E_a;
        float term_2=m_B*(m_b+m_B)*Q+m_B*(m_B+m_b-m_a)*E_a;
        
        float gamma = sqrt(term_1/term_2);
        //cout<<gamma<<endl;
        
        float theta=TMath::ACos((gamma+cos(theta_cm))/sqrt(1+gamma*gamma+2*gamma*cos(theta_cm)));
        
        float terma=(m_B+m_b);
        float termb=-2*sqrt(E_a*m_a*m_b)*cos(theta);
        float termc= - E_a*(m_B-m_a)-Q*m_B;
        
        float sqrtE_b_1=(-termb+sqrt(termb*termb-4*terma*termc))/2/terma;
        float sqrtE_b_2=(-termb-sqrt(termb*termb-4*terma*termc))/2/terma;
        float E_b_1=0;
        float E_b_2=0;
        float E_B_1=0;
        float E_B_2=0;
        
        if(sqrtE_b_1>=0){
            E_b_1=sqrtE_b_1*sqrtE_b_1;
            E_B_1=E_a+Q-E_b_1;
        }
        
        if(sqrtE_b_2>=0){
            E_b_2=sqrtE_b_2*sqrtE_b_2;
            E_B_2=E_a+Q-E_b_2;
        }
        
        cout<<theta_cm/radcon<<" "<<theta/radcon<<" "<<E_b_1<<" "<<E_B_1<<" "<<E_b_2<<" "<<E_B_2<<endl;
    }

}