UserArea/StatTools/qStar.cc

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

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"
#include "xsecLimits.C"

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

  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;

  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;

  pydat3_.mdme[0][4070] = 1;
  pydat3_.mdme[0][4071] = 1;
  pydat3_.mdme[0][4072] = 1;
  pydat3_.mdme[0][4073] = 1;
  pydat3_.mdme[0][4074] = 1;
  pydat3_.mdme[0][4075] = 1;
  pydat3_.mdme[0][4076] = 1;
  pydat3_.mdme[0][4077] = 1;

//  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<<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<<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: "<<pyint5_.xsec[2][0]<<" and the limit: "<<limitXsec<<endl;

  double closestXsec=0, Lambda=0;
  for(double lambda=500; lambda<1500; lambda++){
    // Cross section of the u*->uZ for various Lambdas
    double xsec1 = 
      xsecU * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the u* production
      *uStar2qZ(mqStar,mqStar)/(qStar2qg(mqStar,mqStar,fs) + uStar2qA(mqStar,mqStar) + qStar2qW(mqStar,mqStar) + uStar2qZ(mqStar,mqStar) + qStar2qqq(mqStar,mqStar) + qStar2qll(mqStar,mqStar)) // u*->uZ
      * 0.03366; // Z->mm
    double xsec2 = 
      xsecD * (mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda)*(mqStar/lambda) // total xsec of the d* production
      *dStar2qZ(mqStar,mqStar)/(qStar2qg(mqStar,mqStar,fs) + dStar2qA(mqStar,mqStar) + qStar2qW(mqStar,mqStar) + dStar2qZ(mqStar,mqStar) + qStar2qqq(mqStar,mqStar) + qStar2qll(mqStar,mqStar)) // u*->uZ
      * 0.03366; // Z->mm

    if(lambda>999 && lambda<1001 ) cout<<"lambda="<<lambda<<" xsec1="<<xsec1<<" xsec2="<<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.1
Committed:	Wed Feb 16 15:25:33 2011 UTC (14 years, 2 months ago) by kkotov
Content type:	text/plain
Branch:	MAIN
Log Message:	* empty log message *
#	User	Rev	Content
1	kkotov	1.1	#include<math.h>
2			#include<stdlib.h>
3			#include<iostream>
4
5			using namespace std;
6
7			extern "C" {
8			void pyinit_(char , char, char, double, int, int, int);
9			void pystat_(int *);
10			void pylist_(int *);
11			void pyevnt_(void);
12			void pyckbd_(void);
13			void initpydata_(void);
14			double pyalps_(double *);
15			double pyalem_(double *);
16
17			extern struct {
18			int msel, mselpd, msub[500], kfin[81][2];
19			double ckin[200];
20			} pysubs_;
21
22			extern struct {
23			int mstu[200];
24			double paru[200];
25			int mstj[200];
26			double parj[200];
27			} pydat1_;
28
29			extern struct {
30			int kchg[4][500];
31			double pmas[4][500], parf[2000], vckm[4][4];
32			} pydat2_;
33
34			extern struct {
35			int mdcy[3][500], mdme[2][8000];
36			double brat[8000];
37			int kfdp[5][8000];
38			} pydat3_;
39
40			extern struct {
41			int mstp[200];
42			double parp[200];
43			int msti[200];
44			double pari[200];
45			} pypars_;
46
47			extern struct {
48			int itcm[100];
49			double rtcm[100];
50			} pytcsm_;
51
52			const int pyjets_maxn = 4000;
53			extern struct {
54			int n, npad, k[5][pyjets_maxn];
55			double p[5][pyjets_maxn], v[5][pyjets_maxn];
56			} pyjets_;
57
58			extern struct {
59			int mrpy[6];
60			double rrpy[100];
61			} pydatr_;
62
63			extern struct {
64			int ngenpd, ngen[3][501];
65			double xsec[3][501];
66			} pyint5_;
67
68			}
69
70			#include "branchings.C"
71			#include "xsecLimits.C"
72
73			int main(int argc, char *argv[]){
74			// Initialize defaults:
75			// double lambda = atof(argv[1]);
76			double mqStar = atof(argv[1]); //atof(argv[2]);
77			double f = 1; //atof(argv[3]);
78			double fp = 1; //atof(argv[4]);
79			double fs = atof(argv[2]); //atof(argv[5]);
80
81			double limitXsec = ciLimits(mqStar)*1E-9;
82
83			initpydata_(); pyckbd_();
84			// Random seed:
85			pydatr_.mrpy[0] = 300;
86			// Q* contact interaction cards:
87			pysubs_.msel = 0;
88			pysubs_.msub[146] = 0; // d*
89			pysubs_.msub[147] = 0; // u*
90			pysubs_.msub[166] = 0; // d* ci
91			pysubs_.msub[167] = 1; // u* ci
92			pydat2_.pmas[0][342] = mqStar; // mass if d*
93			pydat2_.pmas[0][343] = mqStar; // mass of u*
94			pytcsm_.rtcm[41] = mqStar; // Lambda = mass for now, difference is accounted later
95			pytcsm_.rtcm[43] = f; // f
96			pytcsm_.rtcm[44] = fp; // fp
97			pytcsm_.rtcm[45] = fs; // fs
98
99			pydat3_.mdme[0][189] = 1;
100			pydat3_.mdme[0][190] = 1;
101			pydat3_.mdme[0][191] = 1;
102			pydat3_.mdme[0][192] = 1;
103			pydat3_.mdme[0][193] = 1;
104			pydat3_.mdme[0][194] = 1;
105			pydat3_.mdme[0][195] = 1;
106			pydat3_.mdme[0][196] = 1;
107			pydat3_.mdme[0][197] = 1;
108			pydat3_.mdme[0][198] = 1;
109			pydat3_.mdme[0][199] = 1;
110			pydat3_.mdme[0][200] = 1;
111			pydat3_.mdme[0][201] = 1;
112			pydat3_.mdme[0][202] = 1;
113			pydat3_.mdme[0][203] = 1;
114			pydat3_.mdme[0][204] = 1;
115			pydat3_.mdme[0][205] = 1;
116			pydat3_.mdme[0][206] = 1;
117			pydat3_.mdme[0][207] = 1;
118			pydat3_.mdme[0][208] = 1;
119
120			pydat3_.mdme[0][173] = 1;
121			pydat3_.mdme[0][174] = 1;
122			pydat3_.mdme[0][175] = 1;
123			pydat3_.mdme[0][176] = 1;
124			pydat3_.mdme[0][177] = 1;
125			pydat3_.mdme[0][178] = 1;
126			pydat3_.mdme[0][181] = 1;
127			pydat3_.mdme[0][182] = 1;
128			pydat3_.mdme[0][183] = 1;
129			pydat3_.mdme[0][184] = 1;
130			pydat3_.mdme[0][185] = 1;
131			pydat3_.mdme[0][186] = 1;
132			pydat3_.mdme[0][187] = 1;
133
134			pydat3_.mdme[0][4070] = 1;
135			pydat3_.mdme[0][4071] = 1;
136			pydat3_.mdme[0][4072] = 1;
137			pydat3_.mdme[0][4073] = 1;
138			pydat3_.mdme[0][4074] = 1;
139			pydat3_.mdme[0][4075] = 1;
140			pydat3_.mdme[0][4076] = 1;
141			pydat3_.mdme[0][4077] = 1;
142
143			// pysubs_.ckin[0] = 60.;
144			// pysubs_.ckin[1] = 120.;
145			// Pythia's D6T UE cards:
146			pydat1_.mstj[10] = 3;
147			pydat1_.mstj[21] = 2;
148			pydat1_.parj[70] = 10;
149			pypars_.mstp[1] = 1;
150			pypars_.mstp[32] = 0;
151			pypars_.mstp[50] = 10042;
152			pypars_.mstp[51] = 2;
153			pypars_.mstp[80] = 1;
154			pypars_.mstp[81] = 4;
155			pydat1_.mstu[20] = 1;
156			pypars_.parp[81] = 1.8387;
157			pypars_.parp[88] = 1960.;
158			pypars_.parp[82] = 0.5;
159			pypars_.parp[83] = 0.4;
160			pypars_.parp[89] = 0.16;
161			pypars_.parp[66] = 2.5;
162			pypars_.parp[84] = 1.0;
163			pypars_.parp[85] = 1.0;
164			pypars_.parp[61] = 1.25;
165			pypars_.parp[63] = 0.2;
166			pypars_.mstp[90] = 1;
167			pypars_.parp[90] = 2.1;
168			pypars_.parp[92] = 15.0;
169
170			// Pythia ProQ20 UE cards:
171			/* pydat1_.mstu[20] = 1;
172			pydat1_.mstj[21] = 2;
173			pydat1_.parj[70] = 10.;
174			pypars_.mstp[1] = 1;
175			pypars_.mstp[32] = 0;
176			pypars_.mstp[50] = 7;
177			pypars_.mstp[51] = 1;
178			pypars_.mstp[80] = 1;
179			pypars_.mstp[81] = 4;
180			pydat1_.parj[0] = 0.073;
181			pydat1_.parj[1] = 0.2;
182			pydat1_.parj[2] = 0.94;
183			pydat1_.parj[3] = 0.032;
184			pydat1_.parj[10] = 0.31;
185			pydat1_.parj[11] = 0.4;
186			pydat1_.parj[12] = 0.54;
187			pydat1_.parj[24] = 0.63;
188			pydat1_.parj[25] = 0.12;
189			pydat1_.mstj[10] = 5;
190			pydat1_.parj[20] = 0.313;
191			pydat1_.parj[40] = 0.49;
192			pydat1_.parj[41] = 1.2;
193			pydat1_.parj[45] = 1.0;
194			pydat1_.parj[46] = 1.0;
195			pypars_.parp[61] = 2.9;
196			pypars_.parp[63] = 0.14;
197			pypars_.parp[66] = 2.65;
198			pypars_.parp[81] = 1.9;
199			pypars_.parp[82] = 0.83;
200			pypars_.parp[83] = 0.6;
201			pypars_.parp[84] = 0.86;
202			pypars_.parp[85] = 0.93;
203			pypars_.parp[88] = 1800.;
204			pypars_.parp[89] = 0.22;
205			pypars_.mstp[90] = 1;
206			pypars_.parp[90] = 2.1;
207			pypars_.parp[92] = 5.0;
208			*/
209			double energy = 7000.;
210			char opt1[4] = "CMS";
211			char opt2[2] = "p";
212			char opt3[2] = "p";
213
214			pyinit_(opt1,opt2,opt3,&energy,3,1,1);
215			for(int event=0; event<1000; event++) pyevnt_();
216			double xsecU = pyint5_.xsec[2][0];
217			cout<<"xsecU="<<xsecU<<endl;
218
219			pysubs_.msub[166] = 1; // d* ci
220			pysubs_.msub[167] = 0; // u* ci
221
222			pyinit_(opt1,opt2,opt3,&energy,3,1,1);
223			for(int event=0; event<1000; event++) pyevnt_();
224			double xsecD = pyint5_.xsec[2][0];
225			cout<<"xsecD="<<xsecD<<endl;
226
227			// cout<<"!!!!! alphaS="<<aS(mqStar*mqStar)<<endl;
228			// cout<<"!!!!! alphaE="<<aE(mqStar*mqStar)<<endl;
229			// cout<<" (alphaS_91.2="<<aS(91.291.2)<<" alphaS_500="<<aS(500.500.)<<")"<<endl;
230			// int stat=1;
231			// pystat_(&stat);
232			// stat=13;
233			// pylist_(&stat);
234			cout<<" Total cross section: "<<pyint5_.xsec[2][0]<<" and the limit: "<<limitXsec<<endl;
235
236			double closestXsec=0, Lambda=0;
237			for(double lambda=500; lambda<1500; lambda++){
238			// Cross section of the u*->uZ for various Lambdas
239			double xsec1 =
240			xsecU * (mqStar/lambda)(mqStar/lambda)(mqStar/lambda)(mqStar/lambda) // total xsec of the u production
241			uStar2qZ(mqStar,mqStar)/(qStar2qg(mqStar,mqStar,fs) + uStar2qA(mqStar,mqStar) + qStar2qW(mqStar,mqStar) + uStar2qZ(mqStar,mqStar) + qStar2qqq(mqStar,mqStar) + qStar2qll(mqStar,mqStar)) // u->uZ
242			* 0.03366; // Z->mm
243			double xsec2 =
244			xsecD * (mqStar/lambda)(mqStar/lambda)(mqStar/lambda)(mqStar/lambda) // total xsec of the d production
245			dStar2qZ(mqStar,mqStar)/(qStar2qg(mqStar,mqStar,fs) + dStar2qA(mqStar,mqStar) + qStar2qW(mqStar,mqStar) + dStar2qZ(mqStar,mqStar) + qStar2qqq(mqStar,mqStar) + qStar2qll(mqStar,mqStar)) // u->uZ
246			* 0.03366; // Z->mm
247
248			if(lambda>999 && lambda<1001 ) cout<<"lambda="<<lambda<<" xsec1="<<xsec1<<" xsec2="<<xsec2<<endl;
249			double xsec = xsec1 + xsec2;
250			if( fabs(xsec-limitXsec)<fabs(closestXsec-limitXsec) ){
251			closestXsec=xsec;
252			Lambda = lambda;
253			}
254			}
255			cout<<" Partial cross section: "<<closestXsec<<endl;
256			cout<<" Lambda: "<<Lambda<<endl;
257
258			return 0;
259			}