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/
// To run: ./qStar 1000 1

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 ciLimits(double mass){ //mass in GeV/c/c
return 0.430487-1000.*0.115741/mass+1000.*1000.*0.106742/mass/mass;
}
double giLimits(double mass){ //same
return 0.780309-1000.*1.47909/mass+1000.*1000.*0.981498/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     = atof(argv[2]); //atof(argv[5]);

  double limitXsec = ciLimits(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]    = 0;  // u*
  pysubs_.msub[166]    = 0;  // d* ci
  pysubs_.msub[167]    = 1;  // 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)"<<endl;

  pysubs_.msub[166]    = 1;  // d* ci
  pysubs_.msub[167]    = 0;  // u* ci

  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)"<<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
  double closestXsec=0, Lambda=0;
  for(double lambda=500; lambda<3000; lambda++){
    // Cross section of the u*->uZ for various Lambdas
    double 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
    double xsec1 = xsecU * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the u* production
      * br1 * 0.03366; // u*->uZ x Z->mm
    double 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
    double xsec2 = xsecD * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the d* production
      * br2 * 0.03366; // d*->dZ x Z->mm

    if(lambda>999 && lambda<1001 ) cout<<"lambda="<<lambda<<" 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: "<<closestXsec<<endl;
  cout<<" Lambda: "<<Lambda<<endl;

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