UserArea/StatTools/qStar.cc

#include<math.h>
#include<stdlib.h>
#include<iostream>

// To compile with Pythia:
// g++ -g -c -Wall initpydata.f
// g++ -g -c -Wall pythia-6.4.22.f
// g++ -g -c -Wall qStar.cc
// g++ -o qStar qStar.o pythia-6.4.22.o initpydata.o -L./lib/ -lLHAPDF -lgfortran -static
// setenv LHAPATH /raid/raid3/cms/kkotov/Pythia/share/lhapdf/PDFsets/

using namespace std;

extern "C" {
  void pyinit_(char *, char*, char*, double*, int, int, int);
  void pystat_(int *);
  void pylist_(int *);
  void pyevnt_(void);
  void pyckbd_(void);
  void initpydata_(void);
  double pyalps_(double *);
  double pyalem_(double *);

  extern struct {
    int msel, mselpd, msub[500], kfin[81][2];
    double ckin[200];
  } pysubs_;

  extern struct {
    int mstu[200];
    double paru[200];
    int mstj[200];
    double parj[200];
  } pydat1_;

  extern struct {
    int kchg[4][500];
    double pmas[4][500], parf[2000], vckm[4][4];
  } pydat2_;

  extern struct {
    int mdcy[3][500], mdme[2][8000];
    double brat[8000];
    int kfdp[5][8000];
  } pydat3_;

  extern struct {
    int mstp[200];
    double parp[200];
    int msti[200];
    double pari[200];
  } pypars_;

  extern struct {
    int itcm[100];
    double rtcm[100];
  } pytcsm_;

  const int pyjets_maxn = 4000;
  extern struct {
    int n, npad, k[5][pyjets_maxn];
    double p[5][pyjets_maxn], v[5][pyjets_maxn];
  } pyjets_;

  extern struct {
        int mrpy[6];
        double rrpy[100];
  } pydatr_;

  extern struct {
        int ngenpd, ngen[3][501];
        double xsec[3][501];
  } pyint5_;

}

//#include "branchings.C"
double aS(double q2){ return pyalps_(&q2); }
double aE(double q2){ return pyalem_(&q2); }
double qStar2qg (double m,double lambda, double fS){ return m*aS(m*m)/3.*fS*fS*(m/lambda)*(m/lambda); } // Eq. (10)
// Now reminding several Standard Model facts:
//  1) Only left-handed doublets of the SU(2) isospin group couple to the weak bosons
//  2) For the _left-handed_ quarks the (T3,Y) parameters are:
//     u(1/2,1/3); d(-1/2,1/3); c(1/2,1/3); b(-1/2,1/3); t(1/2,1/3); s(-1/2,1/3)
const double qU = 2./3., qD = 1./3.; // this is eq. (7) assuming f=f'=1
double uStar2qA (double m,double lambda){ return m*aE(m*m)/4.*qU*qU*(m/lambda)*(m/lambda); } // Eq. (5) and (7) for u*
double dStar2qA (double m,double lambda){ return m*aE(m*m)/4.*qD*qD*(m/lambda)*(m/lambda); } // Eq. (5) and (7) for d*
// Some constatns:
const double sin2_theta_W=0.23, cos2_theta_W=1-sin2_theta_W, fZ_uStar=0.35, fZ_dStar=0.42, mW=80.4, mZ=91.2;
// Widths to the W/Z bosons:
double qStar2qW (double m,double lambda){ return m*aE(m*m)/8./sin2_theta_W * 1/2.  
                                                 *(m/lambda)*(m/lambda)*(1-mW*mW/m/m)*(1-mW*mW/m/m)
                                                 *(2+mW*mW/m/m); }
double uStar2qZ (double m,double lambda){ return m*aE(m*m)/8./cos2_theta_W/sin2_theta_W * fZ_uStar*fZ_uStar
                                                 *(m/lambda)*(m/lambda)*(1-mZ*mZ/m/m)*(1-mZ*mZ/m/m)
                                                 *(2+mZ*mZ/m/m); }
double dStar2qZ (double m,double lambda){ return m*aE(m*m)/8./cos2_theta_W/sin2_theta_W * fZ_dStar*fZ_dStar
                                                 *(m/lambda)*(m/lambda)*(1-mZ*mZ/m/m)*(1-mZ*mZ/m/m)
                                                 *(2+mZ*mZ/m/m); }
// C.I. only decays; assume q* can decay to any of 3 generations of quarks and leptons
double qStar2qll(double m,double lambda){ return m/96./3.1415927*(m/lambda)*(m/lambda)*(m/lambda)*(m/lambda)*1*(1   +1+1+1+1+1); }
double qStar2qqq(double m,double lambda){ return m/96./3.1415927*(m/lambda)*(m/lambda)*(m/lambda)*(m/lambda)*3*(4/3.+1+1+1+1+1); }

//#include "xsecLimits.C"
double ciLimitsOld(double mass){ //mass in GeV/c/c
return 0.430487-1000.*0.115741/mass+1000.*1000.*0.106742/mass/mass;
}
double giLimitsOld(double mass){ //same
return 0.780309-1000.*1.47909/mass+1000.*1000.*0.981498/mass/mass;
}
double ciLimits(double mass){ //mass in GeV/c/c
return 0.293825-1000.*0.0635844/mass+1000.*1000.*0.0736631/mass/mass;
}
double giLimits(double mass){ //same
return 0.845417-1000.*1.45942/mass+1000.*1000.*0.840098/mass/mass;
}

int main(int argc, char *argv[]){
// Initialize defaults:
//  double lambda = atof(argv[1]);
  double mqStar = atof(argv[1]); //atof(argv[2]);
  double f      = 1; //atof(argv[3]);
  double fp     = 1; //atof(argv[4]);
  double fs     = 1; //atof(argv[2]); //atof(argv[5]);

  int mode = 1;
  if( argv[2][0] == 'c' ) mode = 1; else mode = -1;

  double limitXsec = ( mode>0 ? ciLimits(mqStar)*1E-9 : giLimits(mqStar)*1E-9 );
  //double limitXsec = ( mode>0 ? ciLimitsOld(mqStar)*1E-9 : giLimitsOld(mqStar)*1E-9 );

  initpydata_(); pyckbd_();
// Random seed:
  pydatr_.mrpy[0]      = 300;
// Q* contact interaction cards:
  pysubs_.msel         = 0;
  pysubs_.msub[146]    = 0;  // d*
  pysubs_.msub[147]    = (mode<0?1:0);  // u*
  pysubs_.msub[166]    = 0;  // d* ci
  pysubs_.msub[167]    = (mode>0?1:0);  // u* ci
  pydat2_.pmas[0][342] = mqStar;  // mass if d*
  pydat2_.pmas[0][343] = mqStar;  // mass of u*
  pytcsm_.rtcm[41]     = mqStar;  // Lambda = mass for now, difference is accounted later
  pytcsm_.rtcm[43]     = f;  // f
  pytcsm_.rtcm[44]     = fp; // fp
  pytcsm_.rtcm[45]     = fs; // fs

  // W decays
  pydat3_.mdme[0][189] = 1;
  pydat3_.mdme[0][190] = 1;
  pydat3_.mdme[0][191] = 1;
  pydat3_.mdme[0][192] = 1;
  pydat3_.mdme[0][193] = 1;
  pydat3_.mdme[0][194] = 1;
  pydat3_.mdme[0][195] = 1;
  pydat3_.mdme[0][196] = 1;
  pydat3_.mdme[0][197] = 1;
  pydat3_.mdme[0][198] = 1;
  pydat3_.mdme[0][199] = 1;
  pydat3_.mdme[0][200] = 1;
  pydat3_.mdme[0][201] = 1;
  pydat3_.mdme[0][202] = 1;
  pydat3_.mdme[0][203] = 1;
  pydat3_.mdme[0][204] = 1;
  pydat3_.mdme[0][205] = 1;
  pydat3_.mdme[0][206] = 1;
  pydat3_.mdme[0][207] = 1;
  pydat3_.mdme[0][208] = 1;

  // Z decays
  pydat3_.mdme[0][173] = 1;
  pydat3_.mdme[0][174] = 1;
  pydat3_.mdme[0][175] = 1;
  pydat3_.mdme[0][176] = 1;
  pydat3_.mdme[0][177] = 1;
  pydat3_.mdme[0][178] = 1;
  pydat3_.mdme[0][181] = 1;
  pydat3_.mdme[0][182] = 1;
  pydat3_.mdme[0][183] = 1;
  pydat3_.mdme[0][184] = 1;
  pydat3_.mdme[0][185] = 1;
  pydat3_.mdme[0][186] = 1;
  pydat3_.mdme[0][187] = 1;

  // q* decays
  pydat3_.mdme[0][4070] = 1; // d*->dg
  pydat3_.mdme[0][4071] = 1; // d*->d\gamma 
  pydat3_.mdme[0][4072] = 1; // d*->dZ
  pydat3_.mdme[0][4073] = 1; // d*->uW
  pydat3_.mdme[0][4074] = 1; // u*->ug
  pydat3_.mdme[0][4075] = 1; // u*->u\gamma
  pydat3_.mdme[0][4076] = 1; // u*->uZ
  pydat3_.mdme[0][4077] = 1; // u*->dW

//  pysubs_.ckin[0]      = 60.;
//  pysubs_.ckin[1]      = 120.;
// Pythia's D6T UE cards:
  pydat1_.mstj[10]     = 3;
  pydat1_.mstj[21]     = 2;
  pydat1_.parj[70]     = 10;
  pypars_.mstp[1]      = 1;
  pypars_.mstp[32]     = 0;
  pypars_.mstp[50]     = 10042;
  pypars_.mstp[51]     = 2;
  pypars_.mstp[80]     = 1;
  pypars_.mstp[81]     = 4;
  pydat1_.mstu[20]     = 1;
  pypars_.parp[81]     = 1.8387;
  pypars_.parp[88]     = 1960.;
  pypars_.parp[82]     = 0.5;
  pypars_.parp[83]     = 0.4;
  pypars_.parp[89]     = 0.16;
  pypars_.parp[66]     = 2.5;
  pypars_.parp[84]     = 1.0;
  pypars_.parp[85]     = 1.0;
  pypars_.parp[61]     = 1.25;
  pypars_.parp[63]     = 0.2;
  pypars_.mstp[90]     = 1;
  pypars_.parp[90]     = 2.1;
  pypars_.parp[92]     = 15.0;

// Pythia ProQ20 UE cards:
/*  pydat1_.mstu[20]     = 1;
  pydat1_.mstj[21]     = 2;
  pydat1_.parj[70]     = 10.;
  pypars_.mstp[1]      = 1;
  pypars_.mstp[32]     = 0;
  pypars_.mstp[50]     = 7;
  pypars_.mstp[51]     = 1;
  pypars_.mstp[80]     = 1;
  pypars_.mstp[81]     = 4;
  pydat1_.parj[0]      = 0.073;
  pydat1_.parj[1]      = 0.2;
  pydat1_.parj[2]      = 0.94;
  pydat1_.parj[3]      = 0.032;
  pydat1_.parj[10]     = 0.31;
  pydat1_.parj[11]     = 0.4;
  pydat1_.parj[12]     = 0.54;
  pydat1_.parj[24]     = 0.63;
  pydat1_.parj[25]     = 0.12;
  pydat1_.mstj[10]     = 5;
  pydat1_.parj[20]     = 0.313;
  pydat1_.parj[40]     = 0.49;
  pydat1_.parj[41]     = 1.2;
  pydat1_.parj[45]     = 1.0;
  pydat1_.parj[46]     = 1.0;
  pypars_.parp[61]     = 2.9;
  pypars_.parp[63]     = 0.14;
  pypars_.parp[66]     = 2.65;
  pypars_.parp[81]     = 1.9;
  pypars_.parp[82]     = 0.83;
  pypars_.parp[83]     = 0.6;
  pypars_.parp[84]     = 0.86;
  pypars_.parp[85]     = 0.93;
  pypars_.parp[88]     = 1800.;
  pypars_.parp[89]     = 0.22;
  pypars_.mstp[90]     = 1;
  pypars_.parp[90]     = 2.1;
  pypars_.parp[92]     = 5.0;
*/
  double energy = 7000.;
  char opt1[4] = "CMS";
  char opt2[2] = "p";
  char opt3[2] = "p";

  pyinit_(opt1,opt2,opt3,&energy,3,1,1);
  for(int event=0; event<1000; event++) pyevnt_();
  double xsecU = pyint5_.xsec[2][0];
  cout<<"xsecU="<<xsecU<<" (for L=m=1000, f=1 this should be ~1527.26 pb for C.I. and ~105.86 pb for G.I.)"<<endl;

  pysubs_.msub[166]    = (mode>0?1:0);  // d* ci
  pysubs_.msub[167]    = 0;  // u* ci

  pysubs_.msub[146]    = (mode<0?1:0);  // d*
  pysubs_.msub[147]    = 0;  // u*

  pyinit_(opt1,opt2,opt3,&energy,3,1,1);
  for(int event=0; event<1000; event++) pyevnt_();
  double xsecD = pyint5_.xsec[2][0];
  cout<<"xsecD="<<xsecD<<" (for L=m=1000, f=1 this should be ~512.004 pb for C.I. and ~46.34 pb for G.I.)"<<endl;

//  cout<<"!!!!! alphaS="<<aS(mqStar*mqStar)<<endl;
//  cout<<"!!!!! alphaE="<<aE(mqStar*mqStar)<<endl;
//  cout<<"  (alphaS_91.2="<<aS(91.2*91.2)<<" alphaS_500="<<aS(500.*500.)<<")"<<endl;
//  int stat=1;
//  pystat_(&stat);
//  stat=13;
//  pylist_(&stat);
  cout<<" Total cross section: "<<(xsecU+xsecD)<<" and the limit: "<<limitXsec<<endl;

  // The code below works strictly for the contact interaction
for(fs=0.0; fs<1.01; fs+=0.1){
  double closestXsec=0, Lambda=0;
  for(double lambda=100; lambda<3000; lambda++){
    // Cross section of the u*->uZ for various Lambdas
    double br1 = 0, xsec1 = 0;
    if( mode>0 ){
       br1 = uStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + uStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + uStar2qZ(mqStar,lambda) + qStar2qqq(mqStar,lambda) + qStar2qll(mqStar,lambda)); // u*->uZ
       xsec1 = xsecU * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the u* production
      * br1 * 0.03366; // u*->uZ x Z->mm
    } else {
       br1 = uStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + uStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + uStar2qZ(mqStar,lambda)); // u*->uZ
       xsec1 = xsecU * (mqStar/lambda)*(mqStar/lambda) * fs*fs // total xsec of the u* production
      * br1 * 0.03366; // u*->uZ x Z->mm
    }

    double br2 = 0, xsec2 = 0;
    if( mode>0 ){
       br2 = dStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + dStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + dStar2qZ(mqStar,lambda) + qStar2qqq(mqStar,lambda) + qStar2qll(mqStar,lambda)); // u*->uZ
       xsec2 = xsecD * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the d* production
      * br2 * 0.03366; // d*->dZ x Z->mm
    } else {
       br2 = dStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + dStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + dStar2qZ(mqStar,lambda)); // u*->uZ
       xsec2 = xsecD * (mqStar/lambda)*(mqStar/lambda) * fs*fs // total xsec of the d* production
      * br2 * 0.03366; // d*->dZ x Z->mm
    }

    if(lambda>500 && lambda<501 ) cout<<"lambda="<<lambda<<" fs="<<fs<<" xsec1(br1="<<br1<<")="<<xsec1<<" xsec2(br2="<<br2<<")="<<xsec2<<endl;
    double xsec = xsec1 + xsec2;
    if( fabs(xsec-limitXsec)<fabs(closestXsec-limitXsec) ){
         closestXsec=xsec;
         Lambda = lambda;
    }
  }
  cout<<" Partial cross section(fs="<<fs<<",m="<<mqStar<<"): "<<closestXsec<<endl;
  cout<<" Lambda(fs="<<fs<<",m="<<mqStar<<"): "<<Lambda<<endl;
}
  return 0;
}
Revision:	1.3
Committed:	Thu Mar 31 19:46:45 2011 UTC (14 years, 1 month ago) by kkotov
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V2012-H-02, V2012-01-00, V2011-01-01, V2011-01-00, AnnaDimuon, V01-00-01, V01-00-00, HEAD
Changes since 1.2:	+47 -19 lines
Log Message:	* empty log message *
#	Content
1	#include<math.h>
2	#include<stdlib.h>
3	#include<iostream>
4
5	// To compile with Pythia:
6	// g++ -g -c -Wall initpydata.f
7	// g++ -g -c -Wall pythia-6.4.22.f
8	// g++ -g -c -Wall qStar.cc
9	// g++ -o qStar qStar.o pythia-6.4.22.o initpydata.o -L./lib/ -lLHAPDF -lgfortran -static
10	// setenv LHAPATH /raid/raid3/cms/kkotov/Pythia/share/lhapdf/PDFsets/
11
12	using namespace std;
13
14	extern "C" {
15	void pyinit_(char , char, char, double, int, int, int);
16	void pystat_(int *);
17	void pylist_(int *);
18	void pyevnt_(void);
19	void pyckbd_(void);
20	void initpydata_(void);
21	double pyalps_(double *);
22	double pyalem_(double *);
23
24	extern struct {
25	int msel, mselpd, msub[500], kfin[81][2];
26	double ckin[200];
27	} pysubs_;
28
29	extern struct {
30	int mstu[200];
31	double paru[200];
32	int mstj[200];
33	double parj[200];
34	} pydat1_;
35
36	extern struct {
37	int kchg[4][500];
38	double pmas[4][500], parf[2000], vckm[4][4];
39	} pydat2_;
40
41	extern struct {
42	int mdcy[3][500], mdme[2][8000];
43	double brat[8000];
44	int kfdp[5][8000];
45	} pydat3_;
46
47	extern struct {
48	int mstp[200];
49	double parp[200];
50	int msti[200];
51	double pari[200];
52	} pypars_;
53
54	extern struct {
55	int itcm[100];
56	double rtcm[100];
57	} pytcsm_;
58
59	const int pyjets_maxn = 4000;
60	extern struct {
61	int n, npad, k[5][pyjets_maxn];
62	double p[5][pyjets_maxn], v[5][pyjets_maxn];
63	} pyjets_;
64
65	extern struct {
66	int mrpy[6];
67	double rrpy[100];
68	} pydatr_;
69
70	extern struct {
71	int ngenpd, ngen[3][501];
72	double xsec[3][501];
73	} pyint5_;
74
75	}
76
77	//#include "branchings.C"
78	double aS(double q2){ return pyalps_(&q2); }
79	double aE(double q2){ return pyalem_(&q2); }
80	double qStar2qg (double m,double lambda, double fS){ return maS(mm)/3.fSfS(m/lambda)(m/lambda); } // Eq. (10)
81	// Now reminding several Standard Model facts:
82	// 1) Only left-handed doublets of the SU(2) isospin group couple to the weak bosons
83	// 2) For the _left-handed_ quarks the (T3,Y) parameters are:
84	// u(1/2,1/3); d(-1/2,1/3); c(1/2,1/3); b(-1/2,1/3); t(1/2,1/3); s(-1/2,1/3)
85	const double qU = 2./3., qD = 1./3.; // this is eq. (7) assuming f=f'=1
86	double uStar2qA (double m,double lambda){ return maE(mm)/4.qUqU(m/lambda)(m/lambda); } // Eq. (5) and (7) for u*
87	double dStar2qA (double m,double lambda){ return maE(mm)/4.qDqD(m/lambda)(m/lambda); } // Eq. (5) and (7) for d*
88	// Some constatns:
89	const double sin2_theta_W=0.23, cos2_theta_W=1-sin2_theta_W, fZ_uStar=0.35, fZ_dStar=0.42, mW=80.4, mZ=91.2;
90	// Widths to the W/Z bosons:
91	double qStar2qW (double m,double lambda){ return maE(mm)/8./sin2_theta_W * 1/2.
92	(m/lambda)(m/lambda)(1-mWmW/m/m)(1-mWmW/m/m)
93	(2+mWmW/m/m); }
94	double uStar2qZ (double m,double lambda){ return maE(mm)/8./cos2_theta_W/sin2_theta_W * fZ_uStar*fZ_uStar
95	(m/lambda)(m/lambda)(1-mZmZ/m/m)(1-mZmZ/m/m)
96	(2+mZmZ/m/m); }
97	double dStar2qZ (double m,double lambda){ return maE(mm)/8./cos2_theta_W/sin2_theta_W * fZ_dStar*fZ_dStar
98	(m/lambda)(m/lambda)(1-mZmZ/m/m)(1-mZmZ/m/m)
99	(2+mZmZ/m/m); }
100	// C.I. only decays; assume q* can decay to any of 3 generations of quarks and leptons
101	double qStar2qll(double m,double lambda){ return m/96./3.1415927(m/lambda)(m/lambda)(m/lambda)(m/lambda)1(1 +1+1+1+1+1); }
102	double qStar2qqq(double m,double lambda){ return m/96./3.1415927(m/lambda)(m/lambda)(m/lambda)(m/lambda)3(4/3.+1+1+1+1+1); }
103
104	//#include "xsecLimits.C"
105	double ciLimitsOld(double mass){ //mass in GeV/c/c
106	return 0.430487-1000.0.115741/mass+1000.1000.*0.106742/mass/mass;
107	}
108	double giLimitsOld(double mass){ //same
109	return 0.780309-1000.1.47909/mass+1000.1000.*0.981498/mass/mass;
110	}
111	double ciLimits(double mass){ //mass in GeV/c/c
112	return 0.293825-1000.0.0635844/mass+1000.1000.*0.0736631/mass/mass;
113	}
114	double giLimits(double mass){ //same
115	return 0.845417-1000.1.45942/mass+1000.1000.*0.840098/mass/mass;
116	}
117
118	int main(int argc, char *argv[]){
119	// Initialize defaults:
120	// double lambda = atof(argv[1]);
121	double mqStar = atof(argv[1]); //atof(argv[2]);
122	double f = 1; //atof(argv[3]);
123	double fp = 1; //atof(argv[4]);
124	double fs = 1; //atof(argv[2]); //atof(argv[5]);
125
126	int mode = 1;
127	if( argv[2][0] == 'c' ) mode = 1; else mode = -1;
128
129	double limitXsec = ( mode>0 ? ciLimits(mqStar)1E-9 : giLimits(mqStar)1E-9 );
130	//double limitXsec = ( mode>0 ? ciLimitsOld(mqStar)1E-9 : giLimitsOld(mqStar)1E-9 );
131
132	initpydata_(); pyckbd_();
133	// Random seed:
134	pydatr_.mrpy[0] = 300;
135	// Q* contact interaction cards:
136	pysubs_.msel = 0;
137	pysubs_.msub[146] = 0; // d*
138	pysubs_.msub[147] = (mode<0?1:0); // u*
139	pysubs_.msub[166] = 0; // d* ci
140	pysubs_.msub[167] = (mode>0?1:0); // u* ci
141	pydat2_.pmas[0][342] = mqStar; // mass if d*
142	pydat2_.pmas[0][343] = mqStar; // mass of u*
143	pytcsm_.rtcm[41] = mqStar; // Lambda = mass for now, difference is accounted later
144	pytcsm_.rtcm[43] = f; // f
145	pytcsm_.rtcm[44] = fp; // fp
146	pytcsm_.rtcm[45] = fs; // fs
147
148	// W decays
149	pydat3_.mdme[0][189] = 1;
150	pydat3_.mdme[0][190] = 1;
151	pydat3_.mdme[0][191] = 1;
152	pydat3_.mdme[0][192] = 1;
153	pydat3_.mdme[0][193] = 1;
154	pydat3_.mdme[0][194] = 1;
155	pydat3_.mdme[0][195] = 1;
156	pydat3_.mdme[0][196] = 1;
157	pydat3_.mdme[0][197] = 1;
158	pydat3_.mdme[0][198] = 1;
159	pydat3_.mdme[0][199] = 1;
160	pydat3_.mdme[0][200] = 1;
161	pydat3_.mdme[0][201] = 1;
162	pydat3_.mdme[0][202] = 1;
163	pydat3_.mdme[0][203] = 1;
164	pydat3_.mdme[0][204] = 1;
165	pydat3_.mdme[0][205] = 1;
166	pydat3_.mdme[0][206] = 1;
167	pydat3_.mdme[0][207] = 1;
168	pydat3_.mdme[0][208] = 1;
169
170	// Z decays
171	pydat3_.mdme[0][173] = 1;
172	pydat3_.mdme[0][174] = 1;
173	pydat3_.mdme[0][175] = 1;
174	pydat3_.mdme[0][176] = 1;
175	pydat3_.mdme[0][177] = 1;
176	pydat3_.mdme[0][178] = 1;
177	pydat3_.mdme[0][181] = 1;
178	pydat3_.mdme[0][182] = 1;
179	pydat3_.mdme[0][183] = 1;
180	pydat3_.mdme[0][184] = 1;
181	pydat3_.mdme[0][185] = 1;
182	pydat3_.mdme[0][186] = 1;
183	pydat3_.mdme[0][187] = 1;
184
185	// q* decays
186	pydat3_.mdme[0][4070] = 1; // d*->dg
187	pydat3_.mdme[0][4071] = 1; // d*->d\gamma
188	pydat3_.mdme[0][4072] = 1; // d*->dZ
189	pydat3_.mdme[0][4073] = 1; // d*->uW
190	pydat3_.mdme[0][4074] = 1; // u*->ug
191	pydat3_.mdme[0][4075] = 1; // u*->u\gamma
192	pydat3_.mdme[0][4076] = 1; // u*->uZ
193	pydat3_.mdme[0][4077] = 1; // u*->dW
194
195	// pysubs_.ckin[0] = 60.;
196	// pysubs_.ckin[1] = 120.;
197	// Pythia's D6T UE cards:
198	pydat1_.mstj[10] = 3;
199	pydat1_.mstj[21] = 2;
200	pydat1_.parj[70] = 10;
201	pypars_.mstp[1] = 1;
202	pypars_.mstp[32] = 0;
203	pypars_.mstp[50] = 10042;
204	pypars_.mstp[51] = 2;
205	pypars_.mstp[80] = 1;
206	pypars_.mstp[81] = 4;
207	pydat1_.mstu[20] = 1;
208	pypars_.parp[81] = 1.8387;
209	pypars_.parp[88] = 1960.;
210	pypars_.parp[82] = 0.5;
211	pypars_.parp[83] = 0.4;
212	pypars_.parp[89] = 0.16;
213	pypars_.parp[66] = 2.5;
214	pypars_.parp[84] = 1.0;
215	pypars_.parp[85] = 1.0;
216	pypars_.parp[61] = 1.25;
217	pypars_.parp[63] = 0.2;
218	pypars_.mstp[90] = 1;
219	pypars_.parp[90] = 2.1;
220	pypars_.parp[92] = 15.0;
221
222	// Pythia ProQ20 UE cards:
223	/* pydat1_.mstu[20] = 1;
224	pydat1_.mstj[21] = 2;
225	pydat1_.parj[70] = 10.;
226	pypars_.mstp[1] = 1;
227	pypars_.mstp[32] = 0;
228	pypars_.mstp[50] = 7;
229	pypars_.mstp[51] = 1;
230	pypars_.mstp[80] = 1;
231	pypars_.mstp[81] = 4;
232	pydat1_.parj[0] = 0.073;
233	pydat1_.parj[1] = 0.2;
234	pydat1_.parj[2] = 0.94;
235	pydat1_.parj[3] = 0.032;
236	pydat1_.parj[10] = 0.31;
237	pydat1_.parj[11] = 0.4;
238	pydat1_.parj[12] = 0.54;
239	pydat1_.parj[24] = 0.63;
240	pydat1_.parj[25] = 0.12;
241	pydat1_.mstj[10] = 5;
242	pydat1_.parj[20] = 0.313;
243	pydat1_.parj[40] = 0.49;
244	pydat1_.parj[41] = 1.2;
245	pydat1_.parj[45] = 1.0;
246	pydat1_.parj[46] = 1.0;
247	pypars_.parp[61] = 2.9;
248	pypars_.parp[63] = 0.14;
249	pypars_.parp[66] = 2.65;
250	pypars_.parp[81] = 1.9;
251	pypars_.parp[82] = 0.83;
252	pypars_.parp[83] = 0.6;
253	pypars_.parp[84] = 0.86;
254	pypars_.parp[85] = 0.93;
255	pypars_.parp[88] = 1800.;
256	pypars_.parp[89] = 0.22;
257	pypars_.mstp[90] = 1;
258	pypars_.parp[90] = 2.1;
259	pypars_.parp[92] = 5.0;
260	*/
261	double energy = 7000.;
262	char opt1[4] = "CMS";
263	char opt2[2] = "p";
264	char opt3[2] = "p";
265
266	pyinit_(opt1,opt2,opt3,&energy,3,1,1);
267	for(int event=0; event<1000; event++) pyevnt_();
268	double xsecU = pyint5_.xsec[2][0];
269	cout<<"xsecU="<<xsecU<<" (for L=m=1000, f=1 this should be ~1527.26 pb for C.I. and ~105.86 pb for G.I.)"<<endl;
270
271	pysubs_.msub[166] = (mode>0?1:0); // d* ci
272	pysubs_.msub[167] = 0; // u* ci
273
274	pysubs_.msub[146] = (mode<0?1:0); // d*
275	pysubs_.msub[147] = 0; // u*
276
277	pyinit_(opt1,opt2,opt3,&energy,3,1,1);
278	for(int event=0; event<1000; event++) pyevnt_();
279	double xsecD = pyint5_.xsec[2][0];
280	cout<<"xsecD="<<xsecD<<" (for L=m=1000, f=1 this should be ~512.004 pb for C.I. and ~46.34 pb for G.I.)"<<endl;
281
282	// cout<<"!!!!! alphaS="<<aS(mqStar*mqStar)<<endl;
283	// cout<<"!!!!! alphaE="<<aE(mqStar*mqStar)<<endl;
284	// cout<<" (alphaS_91.2="<<aS(91.291.2)<<" alphaS_500="<<aS(500.500.)<<")"<<endl;
285	// int stat=1;
286	// pystat_(&stat);
287	// stat=13;
288	// pylist_(&stat);
289	cout<<" Total cross section: "<<(xsecU+xsecD)<<" and the limit: "<<limitXsec<<endl;
290
291	// The code below works strictly for the contact interaction
292	for(fs=0.0; fs<1.01; fs+=0.1){
293	double closestXsec=0, Lambda=0;
294	for(double lambda=100; lambda<3000; lambda++){
295	// Cross section of the u*->uZ for various Lambdas
296	double br1 = 0, xsec1 = 0;
297	if( mode>0 ){
298	br1 = uStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + uStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + uStar2qZ(mqStar,lambda) + qStar2qqq(mqStar,lambda) + qStar2qll(mqStar,lambda)); // u*->uZ
299	xsec1 = xsecU * (mqStar/lambda)(mqStar/lambda)(mqStar/lambda)(mqStar/lambda) // total xsec of the u production
300	* br1 * 0.03366; // u*->uZ x Z->mm
301	} else {
302	br1 = uStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + uStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + uStar2qZ(mqStar,lambda)); // u*->uZ
303	xsec1 = xsecU * (mqStar/lambda)(mqStar/lambda) fsfs // total xsec of the u production
304	* br1 * 0.03366; // u*->uZ x Z->mm
305	}
306
307	double br2 = 0, xsec2 = 0;
308	if( mode>0 ){
309	br2 = dStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + dStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + dStar2qZ(mqStar,lambda) + qStar2qqq(mqStar,lambda) + qStar2qll(mqStar,lambda)); // u*->uZ
310	xsec2 = xsecD * (mqStar/lambda)(mqStar/lambda)(mqStar/lambda)(mqStar/lambda) // total xsec of the d production
311	* br2 * 0.03366; // d*->dZ x Z->mm
312	} else {
313	br2 = dStar2qZ(mqStar,lambda)/(qStar2qg(mqStar,lambda,fs) + dStar2qA(mqStar,lambda) + qStar2qW(mqStar,lambda) + dStar2qZ(mqStar,lambda)); // u*->uZ
314	xsec2 = xsecD * (mqStar/lambda)(mqStar/lambda) fsfs // total xsec of the d production
315	* br2 * 0.03366; // d*->dZ x Z->mm
316	}
317
318	if(lambda>500 && lambda<501 ) cout<<"lambda="<<lambda<<" fs="<<fs<<" xsec1(br1="<<br1<<")="<<xsec1<<" xsec2(br2="<<br2<<")="<<xsec2<<endl;
319	double xsec = xsec1 + xsec2;
320	if( fabs(xsec-limitXsec)<fabs(closestXsec-limitXsec) ){
321	closestXsec=xsec;
322	Lambda = lambda;
323	}
324	}
325	cout<<" Partial cross section(fs="<<fs<<",m="<<mqStar<<"): "<<closestXsec<<endl;
326	cout<<" Lambda(fs="<<fs<<",m="<<mqStar<<"): "<<Lambda<<endl;
327	}
328	return 0;
329	}