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


class drellyan: public jets {

  /* computes the rate for ppbar -> mu+mu- from 
                  parton+parton collisions    */

 public:

  drellyan(double ecm, double ptmin): jets(ecm,ptmin), 
                      P("mrst/lo2002.dat"){}

  mrst P;

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

  double qfcrosssection(double shat, double that, double uhat, int f){
    double Qf = - 1.0/3.0;
    double I3f = -0.5;
    if (f == 2) {
       Qf = 2.0/3.0; 
       I3f = 0.5;
    }
    double alpha = 1.0/129.0;
    double sstw = 0.23;
    double cstw = 0.77;
    double MZ = 91.188;
    double GZ = 2.50; 
    complex Zfactor = shat/(shat - MZ *(MZ - I*GZ));
    complex GLL =  Qf  + Zfactor * (0.5 - sstw)*(I3f - sstw*Qf)/(sstw*cstw);
    complex GLR =  Qf  - Zfactor * (I3f - sstw*Qf)/cstw; 
    complex GRL =  Qf  - Zfactor * (0.5 - sstw)* Qf/cstw;
    complex GRR =  Qf  + Zfactor * sstw * Qf/cstw;
    return  (PI* alpha * alpha * pb /(6.0 *shat))   
                *( (uhat*uhat/(shat*shat))* ( norm(GLL) + norm(GRR)) 
                   +  (that*that/(shat*shat))* ( norm(GLR) + norm(GRL)) );
  }

  double crosssection(){
    if (shat < 50.0) return 0.0;    
      /* with this statement, we can take a small min pt */
                    /* parton densities  */
    double fu1 = P.xfu(x1,Q)/x1;
    double fd1 = P.xfd(x1,Q)/x1;
    double fs1 = P.xfs(x1,Q)/x1;
    double fc1 = P.xfc(x1,Q)/x1;
    double fb1 = P.xfb(x1,Q)/x1;
    double fub1 = P.xfubar(x1,Q)/x1;
    double fdb1 = P.xfdbar(x1,Q)/x1;
    double fsb1 = P.xfsbar(x1,Q)/x1;
    double fcb1 = P.xfcbar(x1,Q)/x1;
    double fbb1 = P.xfbbar(x1,Q)/x1;
    double fu2 = P.xfu(x2,Q)/x2;
    double fd2 = P.xfd(x2,Q)/x2;
    double fs2 = P.xfs(x2,Q)/x2;
    double fc2 = P.xfc(x2,Q)/x2;
    double fb2 = P.xfb(x2,Q)/x2;
    double fub2 = P.xfubar(x2,Q)/x2;
    double fdb2 = P.xfdbar(x2,Q)/x2;
    double fsb2 = P.xfsbar(x2,Q)/x2;
    double fcb2 = P.xfcbar(x2,Q)/x2;
    double fbb2 = P.xfbbar(x2,Q)/x2;
    double fg1 = P.xfg(x1,Q)/x1;
    double fg2 = P.xfg(x2,Q)/x2;
    double qqcs = (fd1*fd2 + fs1*fs2 + fb1*fb2
                + fdb1*fdb2 + fsb1*fsb2 + fbb1*fbb2)
                            *  qfcrosssection(shat,that,uhat,1) 
                  + (fu1*fu2 + fc1*fc2 + fub1*fub2 + fcb1*fcb2)
		            *  qfcrosssection(shat,that,uhat,2);
     return qqcs;
  }

};


main(){

   gnuplot G1("dyptplot.gp",20.0, 150.0, 10.0, 1000000.0);
   G1.setylog();
   G1.bottom(" pt (GeV) ");
   G1.left( " d sigma/ d pt   (pb/GeV) ");
   G1.begin();

   gnuplot G2("dymjjplot.gp",50.0, 400.0, 10.0, 1000000.0);
   G2.setylog();
   G2.bottom(" M mu+mu- (GeV) ");
   G2.left( " d sigma/ d M(mu+mu-)   (pb/GeV) ");
   G2.begin();

   drellyan DY(2000.0, 10.0);  // small min pt; but we constrain shat above

   DY.prepare(1000000);

   double binsize = 2.0;
   histogram H1(10.0, 400.0, binsize);    
   histogram H2(40.0, 1000.0, binsize);    
        /* histogram from 50 to 1000 in bins of 2 GeV */
  
 
     /* sample the distribution of produced jets  */
 
   int nevents = 2000000;
   double weight, error; 

   for (int n=1; n <= nevents; n++){
       if (n%100000 == 0) cout << " getting event " << n << endl;
       /* generate a sample point; record in the histogram with 
              the correct weight   */
       DVector X = DY.getPoint(weight);
       double Mjj = 2.0 * DY.p;
       double pt =  DY.pt;
       H1.add(pt,weight);   
       H2.add(Mjj,weight);   
   }

   double totalcs = DY.integral(error);
       /* total cross section */

   G1.printScaledHistogram(H1,totalcs/binsize,3);
   G2.printScaledHistogram(H2,totalcs/binsize,3);

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