UserArea/StatTools/qStarCI.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);

  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_;
}

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

  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]    = 1;  // u* ci
  pysubs_.msub[167]    = 1;  // d* ci
  pydat2_.pmas[0][342] = mqStar;  // mass if d*
  pydat2_.pmas[0][343] = mqStar;  // mass of u*
  pytcsm_.rtcm[41]     = lambda;  // Lambda = mass
  pytcsm_.rtcm[43]     = f;  // f
  pytcsm_.rtcm[44]     = fp; // fp
  pytcsm_.rtcm[45]     = fs; // fs
  pydat3_.mdme[0][288] = 0;
  pydat3_.mdme[0][289] = 0;
  pydat3_.mdme[0][290] = 0;
  pydat3_.mdme[0][291] = 0;
  pydat3_.mdme[0][292] = 0;
  pydat3_.mdme[0][293] = 0;
  pydat3_.mdme[0][294] = 0;
  pydat3_.mdme[0][295] = 0;
  pydat3_.mdme[0][296] = 0;
  pydat3_.mdme[0][297] = 0;
  pydat3_.mdme[0][298] = 0;
  pydat3_.mdme[0][299] = 0;
  pydat3_.mdme[0][300] = 0;
  pydat3_.mdme[0][301] = 0;
  pydat3_.mdme[0][302] = 0;
  pydat3_.mdme[0][303] = 0;
  pydat3_.mdme[0][304] = 0;
  pydat3_.mdme[0][305] = 0;
  pydat3_.mdme[0][306] = 0;
  pydat3_.mdme[0][307] = 1; // 1
  pydat3_.mdme[0][308] = 0;
  pydat3_.mdme[0][309] = 0;

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

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

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

//  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_();
    if(event%1000==0) cout<<"Event = "<<event<<endl;
/*    int index=0;
    for(index=0; index<pyjets_.n; index++)
      if( pyjets_.k[2][index]==0 && abs(pyjets_.k[1][index])==23 ) break; // Original Z0
    kine.px = pyjets_.p[0][index];
    kine.py = pyjets_.p[1][index];
    kine.pt = sqrt(kine.px*kine.px+kine.py*kine.py);
    kine.pz = pyjets_.p[2][index];
    kine.E  = pyjets_.p[3][index];
    kine.m  = pyjets_.p[4][index];
    dy.Fill();
    pt.Fill(kine.pt);
*/  }
  int stat=1;
  pystat_(&stat);

//  stat=13;
//  pylist_(&stat);
  return 0;
}
Revision:	1.1
Committed:	Wed Dec 8 10:57:12 2010 UTC (14 years, 4 months 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
Log Message:	* empty log message *
#	Content
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
15	extern struct {
16	int msel, mselpd, msub[500], kfin[81][2];
17	double ckin[200];
18	} pysubs_;
19
20	extern struct {
21	int mstu[200];
22	double paru[200];
23	int mstj[200];
24	double parj[200];
25	} pydat1_;
26
27	extern struct {
28	int kchg[4][500];
29	double pmas[4][500], parf[2000], vckm[4][4];
30	} pydat2_;
31
32	extern struct {
33	int mdcy[3][500], mdme[2][8000];
34	double brat[8000];
35	int kfdp[5][8000];
36	} pydat3_;
37
38	extern struct {
39	int mstp[200];
40	double parp[200];
41	int msti[200];
42	double pari[200];
43	} pypars_;
44
45	extern struct {
46	int itcm[100];
47	double rtcm[100];
48	} pytcsm_;
49
50	const int pyjets_maxn = 4000;
51	extern struct {
52	int n, npad, k[5][pyjets_maxn];
53	double p[5][pyjets_maxn], v[5][pyjets_maxn];
54	} pyjets_;
55
56	extern struct {
57	int mrpy[6];
58	double rrpy[100];
59	} pydatr_;
60	}
61
62	int main(int argc, char *argv[]){
63	// Initialize defaults:
64	double lambda = atof(argv[1]);
65	double mqStar = atof(argv[2]);
66	double f = atof(argv[3]);
67	double fp = atof(argv[4]);
68	double fs = atof(argv[5]);
69
70	initpydata_(); pyckbd_();
71	// Random seed:
72	pydatr_.mrpy[0] = 300;
73	// Q* contact interaction cards:
74	pysubs_.msel = 0;
75	pysubs_.msub[146] = 0; // d*
76	pysubs_.msub[147] = 0; // u*
77	pysubs_.msub[166] = 1; // u* ci
78	pysubs_.msub[167] = 1; // d* ci
79	pydat2_.pmas[0][342] = mqStar; // mass if d*
80	pydat2_.pmas[0][343] = mqStar; // mass of u*
81	pytcsm_.rtcm[41] = lambda; // Lambda = mass
82	pytcsm_.rtcm[43] = f; // f
83	pytcsm_.rtcm[44] = fp; // fp
84	pytcsm_.rtcm[45] = fs; // fs
85	pydat3_.mdme[0][288] = 0;
86	pydat3_.mdme[0][289] = 0;
87	pydat3_.mdme[0][290] = 0;
88	pydat3_.mdme[0][291] = 0;
89	pydat3_.mdme[0][292] = 0;
90	pydat3_.mdme[0][293] = 0;
91	pydat3_.mdme[0][294] = 0;
92	pydat3_.mdme[0][295] = 0;
93	pydat3_.mdme[0][296] = 0;
94	pydat3_.mdme[0][297] = 0;
95	pydat3_.mdme[0][298] = 0;
96	pydat3_.mdme[0][299] = 0;
97	pydat3_.mdme[0][300] = 0;
98	pydat3_.mdme[0][301] = 0;
99	pydat3_.mdme[0][302] = 0;
100	pydat3_.mdme[0][303] = 0;
101	pydat3_.mdme[0][304] = 0;
102	pydat3_.mdme[0][305] = 0;
103	pydat3_.mdme[0][306] = 0;
104	pydat3_.mdme[0][307] = 1; // 1
105	pydat3_.mdme[0][308] = 0;
106	pydat3_.mdme[0][309] = 0;
107
108	pydat3_.mdme[0][189] = 0;
109	pydat3_.mdme[0][190] = 0;
110	pydat3_.mdme[0][191] = 0;
111	pydat3_.mdme[0][192] = 0;
112	pydat3_.mdme[0][193] = 0;
113	pydat3_.mdme[0][194] = 0;
114	pydat3_.mdme[0][195] = 0;
115	pydat3_.mdme[0][196] = 0;
116	pydat3_.mdme[0][197] = 0;
117	pydat3_.mdme[0][198] = 0;
118	pydat3_.mdme[0][199] = 0;
119	pydat3_.mdme[0][200] = 0;
120	pydat3_.mdme[0][201] = 0;
121	pydat3_.mdme[0][202] = 0;
122	pydat3_.mdme[0][203] = 0;
123	pydat3_.mdme[0][204] = 0;
124	pydat3_.mdme[0][205] = 0;
125	pydat3_.mdme[0][206] = 0;
126	pydat3_.mdme[0][207] = 0;
127	pydat3_.mdme[0][208] = 0;
128
129	pydat3_.mdme[0][173] = 0;
130	pydat3_.mdme[0][174] = 0;
131	pydat3_.mdme[0][175] = 0;
132	pydat3_.mdme[0][176] = 0;
133	pydat3_.mdme[0][177] = 0;
134	pydat3_.mdme[0][178] = 0;
135	pydat3_.mdme[0][181] = 0;
136	pydat3_.mdme[0][182] = 0;
137	pydat3_.mdme[0][183] = 1;
138	pydat3_.mdme[0][184] = 0;
139	pydat3_.mdme[0][185] = 0;
140	pydat3_.mdme[0][186] = 0;
141	pydat3_.mdme[0][187] = 0;
142
143	pydat3_.mdme[0][4070] = 0;
144	pydat3_.mdme[0][4071] = 0;
145	pydat3_.mdme[0][4072] = 1;
146	pydat3_.mdme[0][4073] = 0;
147	pydat3_.mdme[0][4074] = 0;
148	pydat3_.mdme[0][4075] = 0;
149	pydat3_.mdme[0][4076] = 1;
150	pydat3_.mdme[0][4077] = 0;
151
152	// pysubs_.ckin[0] = 60.;
153	// pysubs_.ckin[1] = 120.;
154	// Pythia's D6T UE cards:
155	pydat1_.mstj[10] = 3;
156	pydat1_.mstj[21] = 2;
157	pydat1_.parj[70] = 10;
158	pypars_.mstp[1] = 1;
159	pypars_.mstp[32] = 0;
160	pypars_.mstp[50] = 10042;
161	pypars_.mstp[51] = 2;
162	pypars_.mstp[80] = 1;
163	pypars_.mstp[81] = 4;
164	pydat1_.mstu[20] = 1;
165	pypars_.parp[81] = 1.8387;
166	pypars_.parp[88] = 1960.;
167	pypars_.parp[82] = 0.5;
168	pypars_.parp[83] = 0.4;
169	pypars_.parp[89] = 0.16;
170	pypars_.parp[66] = 2.5;
171	pypars_.parp[84] = 1.0;
172	pypars_.parp[85] = 1.0;
173	pypars_.parp[61] = 1.25;
174	pypars_.parp[63] = 0.2;
175	pypars_.mstp[90] = 1;
176	pypars_.parp[90] = 2.1;
177	pypars_.parp[92] = 15.0;
178
179	// Pythia ProQ20 UE cards:
180	/* pydat1_.mstu[20] = 1;
181	pydat1_.mstj[21] = 2;
182	pydat1_.parj[70] = 10.;
183	pypars_.mstp[1] = 1;
184	pypars_.mstp[32] = 0;
185	pypars_.mstp[50] = 7;
186	pypars_.mstp[51] = 1;
187	pypars_.mstp[80] = 1;
188	pypars_.mstp[81] = 4;
189	pydat1_.parj[0] = 0.073;
190	pydat1_.parj[1] = 0.2;
191	pydat1_.parj[2] = 0.94;
192	pydat1_.parj[3] = 0.032;
193	pydat1_.parj[10] = 0.31;
194	pydat1_.parj[11] = 0.4;
195	pydat1_.parj[12] = 0.54;
196	pydat1_.parj[24] = 0.63;
197	pydat1_.parj[25] = 0.12;
198	pydat1_.mstj[10] = 5;
199	pydat1_.parj[20] = 0.313;
200	pydat1_.parj[40] = 0.49;
201	pydat1_.parj[41] = 1.2;
202	pydat1_.parj[45] = 1.0;
203	pydat1_.parj[46] = 1.0;
204	pypars_.parp[61] = 2.9;
205	pypars_.parp[63] = 0.14;
206	pypars_.parp[66] = 2.65;
207	pypars_.parp[81] = 1.9;
208	pypars_.parp[82] = 0.83;
209	pypars_.parp[83] = 0.6;
210	pypars_.parp[84] = 0.86;
211	pypars_.parp[85] = 0.93;
212	pypars_.parp[88] = 1800.;
213	pypars_.parp[89] = 0.22;
214	pypars_.mstp[90] = 1;
215	pypars_.parp[90] = 2.1;
216	pypars_.parp[92] = 5.0;
217	*/
218	double energy = 7000.;
219	char opt1[4] = "CMS";
220	char opt2[2] = "p";
221	char opt3[2] = "p";
222
223	pyinit_(opt1,opt2,opt3,&energy,3,1,1);
224
225	for(int event=0; event<1000; event++){
226	pyevnt_();
227	if(event%1000==0) cout<<"Event = "<<event<<endl;
228	/* int index=0;
229	for(index=0; index<pyjets_.n; index++)
230	if( pyjets_.k[2][index]==0 && abs(pyjets_.k[1][index])==23 ) break; // Original Z0
231	kine.px = pyjets_.p[0][index];
232	kine.py = pyjets_.p[1][index];
233	kine.pt = sqrt(kine.pxkine.px+kine.pykine.py);
234	kine.pz = pyjets_.p[2][index];
235	kine.E = pyjets_.p[3][index];
236	kine.m = pyjets_.p[4][index];
237	dy.Fill();
238	pt.Fill(kine.pt);
239	*/ }
240	int stat=1;
241	pystat_(&stat);
242
243	// stat=13;
244	// pylist_(&stat);
245	return 0;
246	}