#include "GetRate.h"
#include "DIS_Lite.h"
/*GetRate{{{*/
void GetRate(const string spec, const string Com, double E0, double P0, double T0, 
		double* xb_avg,double* W2_avg, double* Q2_avg, double* xs_he3,double* xs_he3_avg,double* rate_he3, double* xs_h3,double* xs_h3_avg,double* rate_h3, double* xs_d2,double* xs_d2_avg,double* rate_d2, double *rate_window){
	bool IsDebug = false;
	//bool IsDebug = true;
	/*
	   double E0 = 11.0; //GeV
	   double P0; //GeV, Scattered Electrons' Central Momentum
	   double T0; //Degree, Scattered Electrons' Central Angle
	   string spec = "";
	   cout<<" --- Which Spectrometer (HRS or BB) = "; cin >> spec;
	   cout<<" --- Beam Energy (E0, GeV) = "; cin >> E0;
	   cout<<" --- Scattered Electron's Central Momentum (P0, GeV) = "; cin >> P0;
	   cout<<" --- Scattered Electron's Central Angle (Theta0, Deg) = "; cin >> T0;
	   */
	int Nevnt = 100000;
	if(IsDebug){cout<<" ----Total Samples = "; cin >> Nevnt;}
	cout<<endl<<endl;

    DIS* dis = new DIS();
    dis->SetKin(E0, P0, T0*Deg2Rad);
 
	Event* evt = new Event();     
	evt->SetSpectrometer(spec);
	if(spec=="BB"){
		if(Com=="narrow")
			evt->SetTGDpRange(-0.045, 0.045);
		else{
			double pmin = (1.0  - P0)/P0;//BB's minimum momentum set to 1GeV
			double pmax = (11.0 - P0)/P0;//BB's maximum momentum set to 11GeV
			evt->SetTGDpRange(pmin, pmax);
		}	
	}	
	double Phase_Space = evt->GetPhaseSpace() * P0; //sr*GeV //it was deltaTheta * deltaPhi *deltaDp

	/*Define Target Info{{{*/
	double Beam_Current = 0.0;//A

	//Upstream Window:   0.010" x 0.069g/cm-2
	//Downstream Window: 0.011" x 0.075g/cm-2
	const double Window_ArialDensity =  ( 0.069+0.075 );//g/cm-2, Al-2219 Alloy
	double Window_Length = 25.;//cm
	if(spec=="HRS")
		Window_Length = 15.0;//cm

	double Target_Density = 0.0;//g*cm-3
	double Target_Thick = 0.0;//cm
	evt->SetVertexRange(-Target_Thick/2.0, Target_Thick/2.0);

	Beam_Current = 20*1e-6;//25uA
	Target_Thick = Window_Length;
	Target_Density = 0.124/25.00;//gcm-3

	double Window_Lumi= Beam_Current/NQe * Window_ArialDensity*Avogadro/27;//cm-2*s-1
	double D2_Target_Lumi= Beam_Current/NQe * Target_Thick*Target_Density*Avogadro/2;//cm-2*s-1
    double sigma0_d2 = dis->Sigma(2,1);

	Beam_Current = 20*1e-6;//25uA
	Target_Thick = Window_Length;
	Target_Density = 0.081/25.00;//g*cm-3

	double H3_Target_Lumi= Beam_Current/NQe * Target_Thick*Target_Density*Avogadro/3;//cm-2*s-1
    double sigma0_h3 = dis->Sigma(3,1);


	Beam_Current = 25*1e-6;//25uA
	Target_Thick = Window_Length;
	Target_Density = 0.093/25.00;//g*cm-3

	double He3_Target_Lumi= Beam_Current/NQe * Target_Thick*Target_Density*Avogadro/3;//cm-2*s-1
    double sigma0_he3 = dis->Sigma(3,2);

	xs_d2[0] = sigma0_d2;
	xs_h3[0] = sigma0_h3;
	xs_he3[0] = sigma0_he3;
	/*}}}*/

	const double CosA = cos( T0*Deg2Rad );
	const double SinA = sin( T0*Deg2Rad );
	double dp_tg,theta_tg, phi_tg,Theta, Ep, sigma_d2, sigma_h3, sigma_he3,sigma_al;
	double xs_p,xs_n,xs_d;
	double xb,nu, Q2, W;   
	double D2_XS_Sum = 0.0,H3_XS_Sum = 0.0,He3_XS_Sum = 0.0;
	double XS_Sum_Window = 0.0;

	/*Start sampling events{{{*/
    int Count=0;  double xb_sum=0.0, W2_sum=0.0, Q2_sum=0.0;
    //randomly generate events within the spectrometer phase space and caluclate the XS
	for(int i=0; i<Nevnt; i++){
		xs_p = -1e-33; xs_n = -1e-33; xs_d = -1e-33; 

		dp_tg = evt->GetTG_Dp();//GeV
		Ep = P0 * (1.0 + dp_tg);//GeV

		theta_tg = evt->GetTG_Theta(); //rad
		phi_tg = evt->GetTG_Phi(); //rad
		Theta = acos( (CosA-phi_tg*SinA) / sqrt(1.+theta_tg*theta_tg+phi_tg*phi_tg) );//rad

		nu = E0-Ep;
		Q2 = 4.0*E0*Ep*pow(sin(Theta/2.0),2);
		xb = Q2 / (2.0*PROTON_MASS*nu);
		W = sqrt( PROTON_MASS*PROTON_MASS + 2.0*PROTON_MASS*nu - Q2 );

		if(IsDebug)
			printf("^^^^ Ep=%8.5fGeV, Theta=%8.5fDeg, xb = %4.3f\n"
					,Ep,Theta/Deg2Rad, xb);

		sigma_d2 = 0;
		sigma_h3 = 0;
		sigma_he3 = 0;
		sigma_al = 0;
        /*Calculate the cross sections{{{*/
        if(xb<1.0){
            dis->SetKin(E0, Ep, Theta);
            xs_p = dis->Sigma("Proton");
            xs_n = dis->Sigma("Neutron");
            xs_d = dis->Sigma("Deutron");

            if(isnan(xs_p)||xs_p<-1e-20)
                xs_p=1e-66;
            if(isnan(xs_n)||xs_n<-1e-20)
                xs_n=1e-66;
            if(isnan(xs_d)||xs_d<-1e-20)
                xs_d=1e-66;

            if(IsDebug)
                printf("^^^^ XS_P = %10.4e, XS_N =  %10.4e\n", xs_p, xs_n);

            sigma_d2 = xs_d;
            sigma_h3 = xs_d+xs_n;
            sigma_he3 = xs_d+xs_p;
            sigma_al = (xs_p * 13 + xs_n * 14);

            if(IsDebug)
                printf("^^^^ XS_D = %10.4e\n", xs_d);

            //if(W>2.0&&xb<1.00){//Only consider the DIS events with W>2.0 cut, for high-x poins, we may need smaller W cut.
            if((W>2.0&&xb<0.70)||(W>sqrt(3.0)&&xb>0.70)){//Only consider the DIS events with W>2.0 cut, for high-x poins, we may need smaller W cut.
                D2_XS_Sum += sigma_d2 * NBARN_TO_CM2;//cm2/GeV/sr
                H3_XS_Sum += sigma_h3 * NBARN_TO_CM2;//cm2/GeV/sr
                He3_XS_Sum += sigma_he3 * NBARN_TO_CM2;//cm2/GeV/sr
                XS_Sum_Window += sigma_al * NBARN_TO_CM2;//Assuming it is pure Al, Z = 13, N = 14
                xb_sum+=xb;
                W2_sum+=W*W;
                Q2_sum+=Q2;
                Count++;
            }
        }
        /*}}}*/

		if(IsDebug)
			printf("XS_Sum=%10.4e, XS_Sum_Window=%10.4e\n",H3_XS_Sum, XS_Sum_Window);
	}
	/*}}}*/

	/*Calculate rates{{{*/
	double He3_Target_Rate = He3_Target_Lumi * He3_XS_Sum * Phase_Space / Nevnt; //cm-2*s-1 * cm2/GeV/sr * sr*GeV = s-1 = Hz
	double H3_Target_Rate = H3_Target_Lumi * H3_XS_Sum * Phase_Space / Nevnt; //cm-2*s-1 * cm2/GeV/sr * sr*GeV = s-1 = Hz
	double D2_Target_Rate = D2_Target_Lumi * D2_XS_Sum * Phase_Space / Nevnt; //cm-2*s-1 * cm2/GeV/sr * sr*GeV = s-1 = Hz
	double Window_Rate = Window_Lumi * XS_Sum_Window * Phase_Space / Nevnt;  //cm-2*s-1 * cm2/GeV/sr * sr*GeV = s-1 = Hz
    
    xb_avg[0] = xb_sum/Count;
    W2_avg[0] = W2_sum/Count;
    Q2_avg[0] = Q2_sum/Count;
	xs_d2_avg[0] = D2_XS_Sum/Count / NBARN_TO_CM2;
	xs_h3_avg[0] = H3_XS_Sum/Count / NBARN_TO_CM2;
	xs_he3_avg[0] = He3_XS_Sum/Count / NBARN_TO_CM2;
	/*}}}*/

	rate_d2[0] = D2_Target_Rate;
	rate_h3[0] = H3_Target_Rate;
	rate_he3[0] = He3_Target_Rate;
	rate_window[0] = Window_Rate;

	delete evt;
    delete dis;
}
/*}}}GetRate*/

/*main{{{*/
int main(){
        string spec;
        const double E0 = 11.0;
        const double Evt_Bin = 25000.;//25K events per data point
    double P0, T0, xbj, Q2, W2;

    double xs_he3, xs_h3, xs_d2, rate_he3, rate_h3, rate_d2;
    double rate_he3_bb, rate_h3_bb, rate_d2_bb, rate_w, rate_w_bb;
    double xs_h3_avg, xs_he3_avg,xs_d2_avg;
    double hour_d2, hour_d2_bb, hour_he3, hour_he3_bb, hour_h3, hour_h3_bb;
    double xbj_avg, W2_avg, Q2_avg;

    cout <<"---- Spectrometer (BB or HRS) "; cin >> spec;
    cout <<"---- E0 = 11~GeV"<<endl;
    cout <<"---- Ep (GeV) = "; cin >> P0;
    cout <<"---- Theta (Degree) = "; cin >> T0;

    double SinSQ = pow(sin(T0*Deg2Rad/2.), 2);
    double nu = E0 - P0;
    Q2 = 4.0 * E0 * P0 *SinSQ;
    W2 = PROTON_MASS*PROTON_MASS + 2.0*PROTON_MASS*nu - Q2;
    xbj = Q2 / 2.0 / nu / PROTON_MASS;

    GetRate(spec.c_str(), "",  E0, P0, T0,&xbj_avg, &W2_avg, &Q2_avg, &xs_he3, &xs_he3_avg, &rate_he3, &xs_h3,&xs_h3_avg, &rate_h3, &xs_d2,&xs_d2_avg, &rate_d2, &rate_w);
    //printf("--- %s: Rate_He3 = %f,  Rate_H3 = %f\n",  spec.c_str(), rate_he3, rate_h3);

    rate_he3_bb = 0.0;
    rate_h3_bb = 0.0;
    rate_d2_bb = 0.0;
    if(spec=="BB"){
        GetRate(spec.c_str(), "narrow",  E0, P0, T0,&xbj_avg,&W2_avg,&Q2_avg, &xs_he3, &xs_he3_avg, &rate_he3_bb, &xs_h3,&xs_h3_avg, &rate_h3_bb, &xs_d2,&xs_d2_avg, &rate_d2_bb, &rate_w_bb);
    }

    //Acount for 10% loss of using SOS quards
    if(spec=="HRS"){
        rate_he3 *=0.9;
        rate_h3 *=0.9;
        rate_d2 *=0.9;
    }

    if(rate_he3>1e-4)
        hour_he3 =  Evt_Bin/(rate_he3*60.*60.);
    else
        hour_he3 = 0;
    if(rate_h3>1e-4)
        hour_h3 =  Evt_Bin/(rate_h3*60.*60.);
    else
        hour_h3 = 0;
    if(rate_d2>1e-4)
        hour_d2 =  Evt_Bin/(rate_d2*60.*60.);
    else
        hour_d2 = 0;
    //if(hour_d2<1.0&&hour_d2>0.00001)  hour_d2 = 1.0;
    //if(hour_h3<1.0&&hour_h3>0.00001)  hour_h3 = 1.0;
    //if(hour_he3<1.0&&hour_h3>0.00001) hour_he3 = 1.0;

    if(rate_he3_bb>1e-4)
        hour_he3_bb =  25000./(rate_he3_bb*60.*60.);
    else
        hour_he3_bb = 0;
    if(rate_h3_bb>1e-4)
        hour_h3_bb =  25000./(rate_h3_bb*60.*60.);
    else
        hour_h3_bb = 0;
    if(rate_d2_bb>1e-4)
        hour_d2_bb =  25000./(rate_d2_bb*60.*60.);
    else
        hour_d2_bb = 0;
    //if(hour_d2_bb<1.0&&hour_d2_bb>0.00001)  hour_d2_bb = 1.0;
    //if(hour_h3_bb<1.0&&hour_h3_bb>0.00001)  hour_h3_bb = 1.0;
    //if(hour_he3_bb<1.0&&hour_h3_bb>0.00001) hour_he3_bb = 1.0;

        cout<<endl<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
    if(spec=="BB") {
        printf("&BigBite: E0 = 11 GeV, Ep = %4.2f GeV/c, Theta = %4.2f Degree\n",P0, T0); 
        printf("          xbj = %4.2f(Avg=%4.2f), W2 = %4.2f(Avg=%4.2f) GeV2, Q2 = %4.2f(Avg=%4.2f) GeV2:\n", xbj,xbj_avg, W2,W2_avg, Q2,Q2_avg); 
        printf("          (For rates, numbers in the brackets are rates with a momentum cut+/-4.5percent)\n");
        printf("--------------------------------------------------------------\n");
        printf("      He3: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f(%5.2f) Hz, Time = %5.2f(%4.1f) hour\n",xs_he3,xs_he3_avg, rate_he3,rate_he3_bb, hour_he3,hour_he3_bb); 
        printf("       H3: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f(%5.2f) Hz, Time = %5.2f(%4.1f) hour\n",xs_h3, xs_h3_avg,rate_h3,rate_h3_bb, hour_h3,hour_h3_bb); 
        printf("       D2: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f(%5.2f) Hz, Time = %5.2f(%4.1f) hour\n",xs_d2, xs_d2_avg,rate_d2,rate_d2_bb, hour_d2,hour_d2_bb); 
        printf("  Windows: Rate = %5.2f(%4.1f) Hz at 20 uA\n", rate_w, rate_w_bb);
        }
    else{
        printf("&&&  HRS: E0 = 11 GeV, Ep = %4.2f GeV/c, Theta = %4.2f Degree\n", P0, T0 ); 
        printf("          xbj = %4.2f(Avg=%4.2f), W2 = %4.2f(Avg=%4.2f) GeV2, Q2 = %4.2f(Avg=%4.2f) GeV2:\n", xbj,xbj_avg, W2,W2_avg, Q2,Q2_avg); 
        printf("--------------------------------------------------------------\n");
        printf("      He3: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f Hz, Time = %4.1f hour\n",xs_he3,xs_he3_avg, rate_he3,hour_he3); 
        printf("       H3: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f Hz, Time = %4.1f hour\n",xs_h3, xs_h3_avg,  rate_h3, hour_h3); 
        printf("       D2: XS = %6.4f(Avg=%6.4f) nb/sr/GeV, Rate = %5.2f Hz, Time = %4.1f hour\n",xs_d2, xs_d2_avg,  rate_d2, hour_d2); 
        printf("  Windows: Rate = %5.2f Hz at 20 uA\n", rate_w);
    }
        cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl<<endl;
}
/*}}}*/