#include "jets.h"
#include "gnuplot.h"


class CDjets: public jets {

  /* computes the rate for heavy boson production at a ppbar collider  */

 public:

  CDjets(double ecm, int Which, double Mass, double gamma): jets(ecm, 5.0),  
           P("mrst/vnvalf1155.dat"), which(Which), M(Mass), Gamma(gamma){}

  mrst P;
  int which;  /*  set which = 0 for D , -1 for Cbar, 1 for C  */
  double M, Gamma;

   /* at the time crosssection is called, x1, x2, etc. are defined */

   /* crosssection is supposed to return   
               (d sigma/ dx1 dx2 d cos theta_CM ) * (2/shat)) 
    for this program, compute   d sigma /dx1 dx2  (independent of cos theta).
      Then the  integral over cos theta will give a factor 2.  If we multiply 
          by   1/shat, the integral will then be the correctly normalized 
                       cross section   */

  double BW(double mys){  
     return   (M * Gamma/PI)/(SQR(mys - M*M) + SQR(M*Gamma)); 
  }

  double crosssection(){
    if (shat < SQR(500.0)) return 0.0;
    const double alphas = 0.12;
    double fu1 = P.xfu(x1,Q)/x1;
    double fd1 = P.xfd(x1,Q)/x1;
    double fub1 = P.xfubar(x1,Q)/x1;
    double fdb1 = P.xfdbar(x1,Q)/x1;
    double fu2 = P.xfu(x2,Q)/x2;
    double fd2 = P.xfd(x2,Q)/x2;
    double fub2 = P.xfubar(x2,Q)/x2;
    double fdb2 = P.xfdbar(x2,Q)/x2;
    double pref =  4.0 * PI * PI * alphas/3.0;
    double cs = 0.0;
    if (which == 1 ) { 
      cs = (fu1 * fub2 + fd1 * fdb2) 
             * (4.0/3.0)* PI * PI * alphas * BW(shat) * pb;
    }   else if (which == -1) { 
      cs = (fub1 * fu2 + fdb1 * fd2) 
             * (4.0/3.0)* PI * PI * alphas * BW(shat) * pb;
    }   else if (which == 0) { 
      cs = (fu1 * fu2  + fd1 * fd2  +  fub1 * fub2 + fdb1 * fdb2)
           *(2.0/9.0)*PI *PI*alphas*BW(shat) * pb;
    }
    return cs/shat;
  }

};

main(){

   double ECM = 2000.0;

   gnuplot G1("massplot.gp",600.0, 1000.0, 0.0, 0.05);
   G1.bottom(" ma(jj) (GeV) ");
   G1.left( " d sigma/ d m(jj)   (pb/GeV) ");
   G1.begin();

   gnuplot G2("Yplot.gp", - 2.0, 2.0, 0.0, 5.0);
   G2.bottom(" Y ");
   G2.left( " d sigma/ d Y  (pb)");
   G2.begin();

   for (int which = -1; which <= 1; which++){

     double Mass = 700.0;
     double Gamma = 28.0;
     if (which == 0){
       Mass = 857.0;
       Gamma = 34.0;
     }
 
     CDjets  CDJ(ECM, which, Mass, Gamma);

     double Mbinsize = 4.0;
     histogram HM(500.0, 1200.0, Mbinsize);    
     double Ybinsize = 0.03;
     histogram HY(-10.0, 10.0, Ybinsize);  

     CDJ.prepare(1000000);  
  
 
   /* sample the distribution of produced jets  */
 
     int nevents = 1000000;
     double weight, error;     

     for (int n=1; n <= nevents; n++){
       if (n%100000 == 0) cout << " getting event " << n << endl;
       DVector X = CDJ.getPoint(weight);
       double Mjj = sqrt(CDJ.x1 * CDJ.x2) * ECM;
       HM.add(Mjj,weight);
       double Y = 0.5* log(CDJ.x1/CDJ.x2);
       HY.add(Y,weight);
          /* generate a sample point; record in the histogram with 
              the correct weight   */
     }

     double totalcs = CDJ.integral(error);
     cout << endl<<   " for case " << which << "  sigma =  " 
                           << totalcs <<" pm " << error  << endl << endl;
          /* total cross section */

     G1.printScaledHistogram(HM,totalcs/Mbinsize,which+2);
     G2.printScaledHistogram(HY,totalcs/Ybinsize,which+2);

   }

   G1.finish();
   G2.finish();
 
   return 0;
}