StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMCAsymMaker.cxx
1 //*-- Author : Renee Fatemi
2 
3 #include "TFile.h"
4 #include "StSpinPool/StMCAsymMaker/StMCAsymMaker.h"
5 #include "StChain.h"
6 #include "StMcEventMaker/StMcEventMaker.h"
7 #include "StMcEventTypes.hh"
8 #include "StMcEvent.hh"
9 #include "St_DataSet.h"
10 #include "St_DataSetIter.h"
11 #include "tables/St_g2t_event_Table.h"
12 #include "tables/St_particle_Table.h"
13 #include "tables/St_g2t_pythia_Table.h"
14 #include "StMuDSTMaker/COMMON/StMuEvent.h"
15 #include "StMuDSTMaker/COMMON/StMuDst.h"
16 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
17 
18 //#include "StMCAsymEvent.h"
19 #include "StPythiaEvent.h"
20 
21 #include "tables/St_particle_Table.h"
22 
23 ClassImp(StMCAsymMaker)
24 
25 StMCAsymMaker::StMCAsymMaker(const char *name):StMaker(name) {
26  mEvent = new StPythiaEvent();
27 }
28 
29 StMCAsymMaker::~StMCAsymMaker() {
30  delete mEvent;
31 }
32 
33 Int_t StMCAsymMaker::Init() {
34  int iset = 0;
35  dssvini2009a_(&iset);
36  return StMaker::Init();
37 }
38 
39 void StMCAsymMaker::Zero() {
40  pid=-10;
41  hard_p= -10;
42  cos_th= -10;
43  x1= -10;
44  x2= -10;
45  s= 0;
46  t= 0;
47  u= 0;
48 
49  partonic_all=0;
50  Q2=0;
51  pid=0;
52 
53  df1_NLO_GSA=0;
54  df2_NLO_GSA=0;
55  weight_NLO_GSA=0;
56 
57  df1_NLO_GSB=0;
58  df2_NLO_GSB=0;
59  weight_NLO_GSB=0;
60 
61  df1_NLO_GSC=0;
62  df2_NLO_GSC=0;
63  weight_NLO_GSC=0;
64 
65  df1_LO=0;
66  df2_LO=0;
67  f1_LO=0;
68  f2_LO=0;
69  weight_LO=0;
70 
71  df1_NLO=0;
72  df2_NLO=0;
73  f1_NLO=0;
74  f2_NLO=0;
75  weight_NLO=0;
76 
77  df1_NLO_g0=0;
78  df2_NLO_g0=0;
79  weight_NLO_g0=0;
80 
81  df1_NLO_gmax=0;
82  df2_NLO_gmax=0;
83  weight_NLO_gmax=0;
84 
85  df1_NLO_gmin=0;
86  df2_NLO_gmin=0;
87  weight_NLO_gmin=0;
88 
89  df1_NLO_m015=0;
90  df2_NLO_m015=0;
91  weight_NLO_m015=0;
92 
93  df1_NLO_m030=0;
94  df2_NLO_m030=0;
95  weight_NLO_m030=0;
96 
97  df1_NLO_m045=0;
98  df2_NLO_m045=0;
99  weight_NLO_m045=0;
100 
101  df1_NLO_m060=0;
102  df2_NLO_m060=0;
103  weight_NLO_m060=0;
104 
105  df1_NLO_m075=0;
106  df2_NLO_m075=0;
107  weight_NLO_m075=0;
108 
109  df1_NLO_m090=0;
110  df2_NLO_m090=0;
111  weight_NLO_m090=0;
112 
113  df1_NLO_m105=0;
114  df2_NLO_m105=0;
115  weight_NLO_m105=0;
116 
117  df1_NLO_p030=0;
118  df2_NLO_p030=0;
119  weight_NLO_p030=0;
120 
121  df1_NLO_p045=0;
122  df2_NLO_p045=0;
123  weight_NLO_p045=0;
124 
125  df1_NLO_p060=0;
126  df2_NLO_p060=0;
127  weight_NLO_p060=0;
128 
129  df1_NLO_p070=0;
130  df2_NLO_p070=0;
131  weight_NLO_p070=0;
132 
133  df1_NLO_DSSV=0;
134  df2_NLO_DSSV=0;
135  weight_NLO_DSSV=0;
136 
137  df1_NLO_DSSV2009a=0;
138  df2_NLO_DSSV2009a=0;
139  weight_NLO_DSSV2009a=0;
140 
141  df1_NLO_LSS1=0;
142  df2_NLO_LSS1=0;
143  weight_NLO_LSS1=0;
144 
145  df1_NLO_LSS2=0;
146  df2_NLO_LSS2=0;
147  weight_NLO_LSS2=0;
148 
149  df1_NLO_LSS3=0;
150  df2_NLO_LSS3=0;
151  weight_NLO_LSS3=0;
152 
153  df1_NLO_LSS2010_delGpos = 0;
154  df2_NLO_LSS2010_delGpos = 0;
155  weight_NLO_LSS2010_delGpos = 0;
156 
157  df1_NLO_LSS2010_chsign_delG = 0;
158  df2_NLO_LSS2010_chsign_delG = 0;
159  weight_NLO_LSS2010_chsign_delG = 0;
160 
161  df1_NLO_AAC1=0;
162  df2_NLO_AAC1=0;
163  weight_NLO_AAC1=0;
164 
165  df1_NLO_AAC2=0;
166  df2_NLO_AAC2=0;
167  weight_NLO_AAC2=0;
168 
169  df1_NLO_AAC3=0;
170  df2_NLO_AAC3=0;
171  weight_NLO_AAC3=0;
172 
173  df1_NLO_BB1=0;
174  df2_NLO_BB1=0;
175  weight_NLO_BB1=0;
176 
177  df1_NLO_BB2=0;
178  df2_NLO_BB2=0;
179  weight_NLO_BB2=0;
180 
181  df1_NLO_BB2010 = 0;
182  df2_NLO_BB2010 = 0;
183  weight_NLO_BB2010 = 0;
184 
185  df1_NLO_DNS1=0;
186  df2_NLO_DNS1=0;
187  weight_NLO_DNS1=0;
188 
189  df1_NLO_DNS2=0;
190  df2_NLO_DNS2=0;
191  weight_NLO_DNS2=0;
192 
193 }
194 
195 void StMCAsymMaker::Clear(const Option_t* c) {
196  Zero();
197  mEvent->Clear(c);
198  StMaker::Clear(c);
199 }
200 //_____________________________________________________________________________
202 Int_t StMCAsymMaker::Make() {
203  //Get StMcEvent to look at GEANT record
204  mcEvent = (StMcEvent*)GetDataSet("StMcEvent");
205  if (!mcEvent) {
206  LOG_WARN << "No StMcEvent" << endm;
207  return kStWarn;
208  }
209 
210  //GET EVTID FROM MuDst
211  muDstMaker = (StMuDstMaker*)GetMaker("MuDst"); assert(muDstMaker);
212  StMuDst* dst = muDstMaker->muDst(); assert(dst);
213  muEvent = dst->event(); assert(muEvent);
214  StEventInfo &info=muEvent->eventInfo();
215  evtid=info.id();
216 
217  //GET GEANT EVENT
218  TDataSet *Event = GetDataSet("geant"); //Event->ls(3);
219 
220  //GET PYTHIA RECORD from particleTable
221  TDataSetIter geantDstI(Event);
222  particleTabPtr = (St_particle *) geantDstI("particle");
223 
224  // if no particle table in the file just skip it
225  if (particleTabPtr!=0){
226 
227  particle_st* particleTable = particleTabPtr->GetTable();//particleTabPtr->Print();
228 
229  //GET EVTID and SUBPROCESS ID from struct g2t_event
230  Pg2t_event=(St_g2t_event *) geantDstI("g2t_event"); //Pg2t_event->Print();
231  g2t_event_st *g2t_event1=Pg2t_event->GetTable();
232  geantID= g2t_event1->n_event;
233  geantPID= g2t_event1->subprocess_id;
234  pid=geantPID;
235  //TEST that geantID==eventID to ensure .geant and .MuDst file are synchronized
236  assert(evtid==geantID);
237 
238  //GET PARTONIC KINEMATICS from struct Pg2t_pythia
239  Pg2t_pythia=(St_g2t_pythia *) geantDstI("g2t_pythia");// Pg2t_pythia->Print();
240  g2t_pythia_st *g2t_pythia1=Pg2t_pythia->GetTable();
241  s= g2t_pythia1-> mand_s;
242  t= g2t_pythia1-> mand_t;
243  u= g2t_pythia1-> mand_u;
244  hard_p= g2t_pythia1->hard_p;
245  cos_th= g2t_pythia1->cos_th;
246  x1= g2t_pythia1->bjor_1;
247  x2= g2t_pythia1->bjor_2;
248 
249  //GET FLAVOR AFTER INTIAL RADIATION BEFORE and AFTER SCATTERING
250  flavor1=particleTable[4].idhep;
251  flavor2=particleTable[5].idhep;
252  flavor3=particleTable[6].idhep;
253  flavor4=particleTable[7].idhep;
254 
255  //GET SCATTERED PARTON RECORD
256  parton1[0]=particleTable[6].idhep;// particle id
257  parton1[1]=particleTable[6].phep[0];//px
258  parton1[2]=particleTable[6].phep[1];//py
259  parton1[3]=particleTable[6].phep[2];//pz
260  parton1[4]=particleTable[6].phep[3];//E
261  parton1[5]=particleTable[6].phep[4];//m
262  parton1[6]=particleTable[6].isthep;//status
263  parton1[7]=particleTable[6].jmohep[0];//moth1
264  parton1[8]=particleTable[6].jmohep[1];//moth2
265  parton1[9]=particleTable[6].jdahep[0];//daughter1
266  parton1[10]=particleTable[6].jdahep[1];//daughter2
267  parton2[0]=particleTable[7].idhep;// particle id
268  parton2[1]=particleTable[7].phep[0];//px
269  parton2[2]=particleTable[7].phep[1];//py
270  parton2[3]=particleTable[7].phep[2];//pz
271  parton2[4]=particleTable[7].phep[3];//E
272  parton2[5]=particleTable[7].phep[4];//m
273  parton2[6]=particleTable[7].isthep;//status
274  parton2[7]=particleTable[7].jmohep[0];//moth1
275  parton2[8]=particleTable[7].jmohep[1];//moth2
276  parton2[9]=particleTable[7].jdahep[0];//daughter1
277  parton2[10]=particleTable[7].jdahep[1];//daughter2
278 
279  if (0){//PRINT OUT PYTHIA RECORD
280  printf("PID/evtid from McEvent = %d,%d; PID/evtid from Table = %d,%d:\n",pid,evtid,geantPID,geantID);
281  printf("row | id | px | py | pz | E | m | status | moth1 | moth2 | daught1 | daught2 |\n");
282  for (int i=0; i<particleTabPtr->GetNRows();++i) {
283  printf(" %d, %d, %f, %f, %f, %f, %f, %d, %d, %d, %d, %d\n",i,particleTable[i].idhep, particleTable[i].phep[0],
284  particleTable[i].phep[1], particleTable[i].phep[2] , particleTable[i].phep[3], particleTable[i].phep[4], particleTable[i].isthep,
285  particleTable[i].jmohep[0], particleTable[i].jmohep[1], particleTable[i].jdahep[0], particleTable[i].jdahep[1]);}
286  cout << endl;
287  cout << "flavor1: " << flavor1 << " flavor2: " << flavor2 << " flavor3: " << flavor3 << " flavor4: " << flavor4 << endl << endl;
288  }
289 
290  //Get partonic a_LL, polarized/unpolarized pdfs using Q2 = partonic_pT^2
291  partonic_all=getPartonicALL(s,t,u,pid,flavor1,flavor2,flavor3,flavor4);
292  Q2=hard_p*hard_p;
293 
294  //NLO GS SCENARIO A
295  df1_NLO_GSA=get_polPDF_NLO_GSA(flavor1,x1,Q2);
296  df2_NLO_GSA=get_polPDF_NLO_GSA(flavor2,x2,Q2);
297  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
298  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
299  weight_NLO_GSA=(df1_NLO_GSA*df2_NLO_GSA*partonic_all)/(f1_NLO*f2_NLO);
300 
301  //NLO GS SCENARIO B
302  df1_NLO_GSB=get_polPDF_NLO_GSB(flavor1,x1,Q2);
303  df2_NLO_GSB=get_polPDF_NLO_GSB(flavor2,x2,Q2);
304  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
305  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
306  weight_NLO_GSB=(df1_NLO_GSB*df2_NLO_GSB*partonic_all)/(f1_NLO*f2_NLO);
307 
308  //NLO GS SCENARIO C
309  df1_NLO_GSC=get_polPDF_NLO_GSC(flavor1,x1,Q2);
310  df2_NLO_GSC=get_polPDF_NLO_GSC(flavor2,x2,Q2);
311  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
312  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
313  weight_NLO_GSC=(df1_NLO_GSC*df2_NLO_GSC*partonic_all)/(f1_NLO*f2_NLO);
314 
315  //LO GRSV
316  df1_LO=get_polPDF_LO(flavor1,x1,Q2);
317  df2_LO=get_polPDF_LO(flavor2,x2,Q2);
318  f1_LO=get_unpolPDF_LO(flavor1,x1,Q2);
319  f2_LO=get_unpolPDF_LO(flavor2,x2,Q2);
320  weight_LO=(df1_LO*df2_LO*partonic_all)/(f1_LO*f2_LO);
321 
322  //NLO GRSV
323  df1_NLO=get_polPDF_NLO(flavor1,x1,Q2);
324  df2_NLO=get_polPDF_NLO(flavor2,x2,Q2);
325  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
326  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
327  weight_NLO=(df1_NLO*df2_NLO*partonic_all)/(f1_NLO*f2_NLO);
328 
329  //NLO_g0 GRSV
330  df1_NLO_g0=get_polPDF_NLO_g0(flavor1,x1,Q2);
331  df2_NLO_g0=get_polPDF_NLO_g0(flavor2,x2,Q2);
332  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
333  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
334  weight_NLO_g0=(df1_NLO_g0*df2_NLO_g0*partonic_all)/(f1_NLO*f2_NLO);
335 
336  //NLO_gmax GRSV
337  df1_NLO_gmax=get_polPDF_NLO_gmax(flavor1,x1,Q2);
338  df2_NLO_gmax=get_polPDF_NLO_gmax(flavor2,x2,Q2);
339  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
340  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
341  weight_NLO_gmax=(df1_NLO_gmax*df2_NLO_gmax*partonic_all)/(f1_NLO*f2_NLO);
342 
343  //NLO_gmin GRSV
344  df1_NLO_gmin=get_polPDF_NLO_gmin(flavor1,x1,Q2);
345  df2_NLO_gmin=get_polPDF_NLO_gmin(flavor2,x2,Q2);
346  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
347  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
348  weight_NLO_gmin=(df1_NLO_gmin*df2_NLO_gmin*partonic_all)/(f1_NLO*f2_NLO);
349 
350  //NLO_m015 GRSV
351  df1_NLO_m015=get_polPDF_NLO_m015(flavor1,x1,Q2);
352  df2_NLO_m015=get_polPDF_NLO_m015(flavor2,x2,Q2);
353  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
354  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
355  weight_NLO_m015=(df1_NLO_m015*df2_NLO_m015*partonic_all)/(f1_NLO*f2_NLO);
356 
357  //NLO_m030 GRSV
358  df1_NLO_m030=get_polPDF_NLO_m030(flavor1,x1,Q2);
359  df2_NLO_m030=get_polPDF_NLO_m030(flavor2,x2,Q2);
360  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
361  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
362  weight_NLO_m030=(df1_NLO_m030*df2_NLO_m030*partonic_all)/(f1_NLO*f2_NLO);
363 
364  //NLO_m045 GRSV
365  df1_NLO_m045=get_polPDF_NLO_m045(flavor1,x1,Q2);
366  df2_NLO_m045=get_polPDF_NLO_m045(flavor2,x2,Q2);
367  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
368  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
369  weight_NLO_m045=(df1_NLO_m045*df2_NLO_m045*partonic_all)/(f1_NLO*f2_NLO);
370 
371  //NLO_m060 GRSV
372  df1_NLO_m060=get_polPDF_NLO_m060(flavor1,x1,Q2);
373  df2_NLO_m060=get_polPDF_NLO_m060(flavor2,x2,Q2);
374  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
375  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
376  weight_NLO_m060=(df1_NLO_m060*df2_NLO_m060*partonic_all)/(f1_NLO*f2_NLO);
377 
378  //NLO_m075 GRSV
379  df1_NLO_m075=get_polPDF_NLO_m075(flavor1,x1,Q2);
380  df2_NLO_m075=get_polPDF_NLO_m075(flavor2,x2,Q2);
381  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
382  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
383  weight_NLO_m075=(df1_NLO_m075*df2_NLO_m075*partonic_all)/(f1_NLO*f2_NLO);
384 
385  //NLO_m090 GRSV
386  df1_NLO_m090=get_polPDF_NLO_m090(flavor1,x1,Q2);
387  df2_NLO_m090=get_polPDF_NLO_m090(flavor2,x2,Q2);
388  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
389  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
390  weight_NLO_m090=(df1_NLO_m090*df2_NLO_m090*partonic_all)/(f1_NLO*f2_NLO);
391 
392  //NLO_m105 GRSV
393  df1_NLO_m105=get_polPDF_NLO_m105(flavor1,x1,Q2);
394  df2_NLO_m105=get_polPDF_NLO_m105(flavor2,x2,Q2);
395  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
396  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
397  weight_NLO_m105=(df1_NLO_m105*df2_NLO_m105*partonic_all)/(f1_NLO*f2_NLO);
398 
399  //NLO_p030 GRSV
400  df1_NLO_p030=get_polPDF_NLO_p030(flavor1,x1,Q2);
401  df2_NLO_p030=get_polPDF_NLO_p030(flavor2,x2,Q2);
402  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
403  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
404  weight_NLO_p030=(df1_NLO_p030*df2_NLO_p030*partonic_all)/(f1_NLO*f2_NLO);
405 
406  //NLO_p045 GRSV
407  df1_NLO_p045=get_polPDF_NLO_p045(flavor1,x1,Q2);
408  df2_NLO_p045=get_polPDF_NLO_p045(flavor2,x2,Q2);
409  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
410  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
411  weight_NLO_p045=(df1_NLO_p045*df2_NLO_p045*partonic_all)/(f1_NLO*f2_NLO);
412 
413  //NLO_p060 GRSV
414  df1_NLO_p060=get_polPDF_NLO_p060(flavor1,x1,Q2);
415  df2_NLO_p060=get_polPDF_NLO_p060(flavor2,x2,Q2);
416  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
417  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
418  weight_NLO_p060=(df1_NLO_p060*df2_NLO_p060*partonic_all)/(f1_NLO*f2_NLO);
419 
420  //NLO_p070 GRSV
421  df1_NLO_p070=get_polPDF_NLO_p070(flavor1,x1,Q2);
422  df2_NLO_p070=get_polPDF_NLO_p070(flavor2,x2,Q2);
423  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
424  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
425  weight_NLO_p070=(df1_NLO_p070*df2_NLO_p070*partonic_all)/(f1_NLO*f2_NLO);
426 
427  //NLO DSSV
428  df1_NLO_DSSV=get_polPDF_NLO_DSSV(flavor1,x1,Q2);
429  df2_NLO_DSSV=get_polPDF_NLO_DSSV(flavor2,x2,Q2);
430  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
431  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
432  weight_NLO_DSSV=(df1_NLO_DSSV*df2_NLO_DSSV*partonic_all)/(f1_NLO*f2_NLO);
433 
434  //NLO DSSV 2009a
435  df1_NLO_DSSV2009a=get_polPDF_NLO_DSSV2009a(flavor1,x1,Q2);
436  df2_NLO_DSSV2009a=get_polPDF_NLO_DSSV2009a(flavor2,x2,Q2);
437  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
438  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
439  weight_NLO_DSSV2009a=(df1_NLO_DSSV2009a/f1_NLO)*(df2_NLO_DSSV2009a/f2_NLO)*partonic_all;
440 
441  //NLO LSS SCENARIO 1
442  df1_NLO_LSS1=get_polPDF_NLO_LSS1(flavor1,x1,Q2);
443  df2_NLO_LSS1=get_polPDF_NLO_LSS1(flavor2,x2,Q2);
444  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
445  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
446  weight_NLO_LSS1=(df1_NLO_LSS1*df2_NLO_LSS1*partonic_all)/(f1_NLO*f2_NLO);
447 
448  //NLO LSS SCENARIO 2
449  df1_NLO_LSS2=get_polPDF_NLO_LSS2(flavor1,x1,Q2);
450  df2_NLO_LSS2=get_polPDF_NLO_LSS2(flavor2,x2,Q2);
451  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
452  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
453  weight_NLO_LSS2=(df1_NLO_LSS2*df2_NLO_LSS2*partonic_all)/(f1_NLO*f2_NLO);
454 
455  //NLO LSS SCENARIO 3
456  df1_NLO_LSS3=get_polPDF_NLO_LSS3(flavor1,x1,Q2);
457  df2_NLO_LSS3=get_polPDF_NLO_LSS3(flavor2,x2,Q2);
458  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
459  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
460  weight_NLO_LSS3=(df1_NLO_LSS3*df2_NLO_LSS3*partonic_all)/(f1_NLO*f2_NLO);
461 
462  //NLO LSS2010 pos x*DeltaG
463  df1_NLO_LSS2010_delGpos = get_polPDF_NLO_LSS2010_delGpos(flavor1,x1,Q2);
464  df2_NLO_LSS2010_delGpos = get_polPDF_NLO_LSS2010_delGpos(flavor2,x2,Q2);
465  f1_NLO = get_unpolPDF_NLO(flavor1,x1,Q2);
466  f2_NLO = get_unpolPDF_NLO(flavor2,x2,Q2);
467  weight_NLO_LSS2010_delGpos = (df1_NLO_LSS2010_delGpos*df2_NLO_LSS2010_delGpos*partonic_all)/(f1_NLO*f2_NLO);
468 
469  //NLO LSS2010 node x*DeltaG
470  df1_NLO_LSS2010_chsign_delG = get_polPDF_NLO_LSS2010_chsign_delG(flavor1,x1,Q2);
471  df2_NLO_LSS2010_chsign_delG = get_polPDF_NLO_LSS2010_chsign_delG(flavor2,x2,Q2);
472  f1_NLO = get_unpolPDF_NLO(flavor1,x1,Q2);
473  f2_NLO = get_unpolPDF_NLO(flavor2,x2,Q2);
474  weight_NLO_LSS2010_chsign_delG = (df1_NLO_LSS2010_chsign_delG*df2_NLO_LSS2010_chsign_delG*partonic_all)/(f1_NLO*f2_NLO);
475 
476  //NLO AAC SCENARIO 1
477  df1_NLO_AAC1=get_polPDF_NLO_AAC1(flavor1,x1,Q2);
478  df2_NLO_AAC1=get_polPDF_NLO_AAC1(flavor2,x2,Q2);
479  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
480  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
481  weight_NLO_AAC1=(df1_NLO_AAC1*df2_NLO_AAC1*partonic_all)/(f1_NLO*f2_NLO);
482 
483  //NLO AAC SCENARIO 2
484  df1_NLO_AAC2=get_polPDF_NLO_AAC2(flavor1,x1,Q2);
485  df2_NLO_AAC2=get_polPDF_NLO_AAC2(flavor2,x2,Q2);
486  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
487  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
488  weight_NLO_AAC2=(df1_NLO_AAC2*df2_NLO_AAC2*partonic_all)/(f1_NLO*f2_NLO);
489 
490  //NLO AAC SCENARIO 3
491  df1_NLO_AAC3=get_polPDF_NLO_AAC3(flavor1,x1,Q2);
492  df2_NLO_AAC3=get_polPDF_NLO_AAC3(flavor2,x2,Q2);
493  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
494  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
495  weight_NLO_AAC3=(df1_NLO_AAC3*df2_NLO_AAC3*partonic_all)/(f1_NLO*f2_NLO);
496 
497  //NLO BB SCENARIO 1
498  df1_NLO_BB1=get_polPDF_NLO_BB1(flavor1,x1,Q2);
499  df2_NLO_BB1=get_polPDF_NLO_BB1(flavor2,x2,Q2);
500  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
501  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
502  weight_NLO_BB1=(df1_NLO_BB1*df2_NLO_BB1*partonic_all)/(f1_NLO*f2_NLO);
503 
504  //NLO BB SCENARIO 2
505  df1_NLO_BB2=get_polPDF_NLO_BB2(flavor1,x1,Q2);
506  df2_NLO_BB2=get_polPDF_NLO_BB2(flavor2,x2,Q2);
507  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
508  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
509  weight_NLO_BB2=(df1_NLO_BB2*df2_NLO_BB2*partonic_all)/(f1_NLO*f2_NLO);
510 
511  //NLO BB2010
512  df1_NLO_BB2010 = get_polPDF_NLO_BB2010(flavor1,x1,Q2);
513  df2_NLO_BB2010 = get_polPDF_NLO_BB2010(flavor2,x2,Q2);
514  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
515  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
516  weight_NLO_BB2010 = (df1_NLO_BB2010*df2_NLO_BB2010*partonic_all)/(f1_NLO*f2_NLO);
517 
518  //NLO DNS SCENARIO 1
519  df1_NLO_DNS1=get_polPDF_NLO_DNS1(flavor1,x1,Q2);
520  df2_NLO_DNS1=get_polPDF_NLO_DNS1(flavor2,x2,Q2);
521  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
522  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
523  weight_NLO_DNS1=(df1_NLO_DNS1*df2_NLO_DNS1*partonic_all)/(f1_NLO*f2_NLO);
524 
525  //NLO DNS SCENARIO 2
526  df1_NLO_DNS2=get_polPDF_NLO_DNS2(flavor1,x1,Q2);
527  df2_NLO_DNS2=get_polPDF_NLO_DNS2(flavor2,x2,Q2);
528  f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
529  f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
530  weight_NLO_DNS2=(df1_NLO_DNS2*df2_NLO_DNS2*partonic_all)/(f1_NLO*f2_NLO);
531 
532  if (0) {
533  printf("LO: df1_LO=%f, df2_LO=%f, f1_LO=%f, f2_LO=%f, weight_LO=%f\n",df1_LO,df2_LO,f1_LO,f2_LO,weight_LO);
534  printf("NLO: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO,df2_NLO,f1_NLO,f2_NLO,weight_NLO);
535  printf("NLO_gmin: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_gmin,df2_NLO_gmin,f1_NLO,f2_NLO,weight_NLO_gmin);
536  printf("NLO_g0: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_g0,df2_NLO_g0,f1_NLO,f2_NLO,weight_NLO_g0);
537  printf("NLO_gmax: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_gmax,df2_NLO_gmax,f1_NLO,f2_NLO,weight_NLO_gmax);
538 
539  printf("NLO_m015: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m015,df2_NLO_m015,f1_NLO,f2_NLO,weight_NLO_m015);
540  printf("NLO_m030: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m030,df2_NLO_m030,f1_NLO,f2_NLO,weight_NLO_m030);
541  printf("NLO_m045: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m045,df2_NLO_m045,f1_NLO,f2_NLO,weight_NLO_m045);
542  printf("NLO_m060: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m060,df2_NLO_m060,f1_NLO,f2_NLO,weight_NLO_m060);
543  printf("NLO_m075: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m075,df2_NLO_m075,f1_NLO,f2_NLO,weight_NLO_m075);
544  printf("NLO_m090: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m090,df2_NLO_m090,f1_NLO,f2_NLO,weight_NLO_m090);
545  printf("NLO_m105: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_m105,df2_NLO_m105,f1_NLO,f2_NLO,weight_NLO_m105);
546 
547  printf("NLO_p030: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_p030,df2_NLO_p030,f1_NLO,f2_NLO,weight_NLO_p030);
548  printf("NLO_p045: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_p045,df2_NLO_p045,f1_NLO,f2_NLO,weight_NLO_p045);
549  printf("NLO_p060: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_p060,df2_NLO_p060,f1_NLO,f2_NLO,weight_NLO_p060);
550  printf("NLO_p070: df1_NLO=%f, df2_NLO=%f, f1_NLO=%f, f2_NLO=%f, weight_NLO=%f\n",df1_NLO_p070,df2_NLO_p070,f1_NLO,f2_NLO,weight_NLO_p070);
551 
552 
553  printf("DSSV: df1_DSSV=%f, df2_DSSV=%f, f1_DSSV=%f, f2_DSSV=%f, weight_DSSV=%f\n",df1_NLO_DSSV,df2_NLO_DSSV,f1_NLO,f2_NLO,weight_NLO_DSSV);
554  printf("LSS1: df1_LSS1=%f, df2_LSS1=%f, f1_LSS1=%f, f2_LSS1=%f, weight_LSS1=%f\n",df1_NLO_LSS1,df2_NLO_LSS1,f1_NLO,f2_NLO,weight_NLO_LSS1);
555  printf("LSS2: df1_LSS2=%f, df2_LSS2=%f, f1_LSS2=%f, f2_LSS2=%f, weight_LSS2=%f\n",df1_NLO_LSS2,df2_NLO_LSS2,f1_NLO,f2_NLO,weight_NLO_LSS2);
556  printf("LSS3: df1_LSS3=%f, df2_LSS3=%f, f1_LSS3=%f, f2_LSS3=%f, weight_LSS3=%f\n",df1_NLO_LSS3,df2_NLO_LSS3,f1_NLO,f2_NLO,weight_NLO_LSS3);
557  printf("AAC1: df1_AAC1=%f, df2_AAC1=%f, f1_AAC1=%f, f2_AAC1=%f, weight_AAC1=%f\n",df1_NLO_AAC1,df2_NLO_AAC1,f1_NLO,f2_NLO,weight_NLO_AAC1);
558  printf("AAC2: df1_AAC2=%f, df2_AAC2=%f, f1_AAC2=%f, f2_AAC2=%f, weight_AAC2=%f\n",df1_NLO_AAC2,df2_NLO_AAC2,f1_NLO,f2_NLO,weight_NLO_AAC2);
559  printf("AAC3: df1_AAC3=%f, df2_AAC3=%f, f1_AAC3=%f, f2_AAC3=%f, weight_AAC3=%f\n",df1_NLO_AAC3,df2_NLO_AAC3,f1_NLO,f2_NLO,weight_NLO_AAC3);
560 
561  printf("BB1: df1_BB1=%f, df2_BB1=%f, f1_BB1=%f, f2_BB1=%f, weight_BB1=%f\n",df1_NLO_BB1,df2_NLO_BB1,f1_NLO,f2_NLO,weight_NLO_BB1);
562  printf("BB2: df1_BB2=%f, df2_BB2=%f, f1_BB2=%f, f2_BB2=%f, weight_BB2=%f\n",df1_NLO_BB2,df2_NLO_BB2,f1_NLO,f2_NLO,weight_NLO_BB2);
563  printf("DNS1: df1_DNS1=%f, df2_DNS1=%f, f1_DNS1=%f, f2_DNS1=%f, weight_DNS1=%f\n",df1_NLO_DNS1,df2_NLO_DNS1,f1_NLO,f2_NLO,weight_NLO_DNS1);
564  printf("DNS2: df1_DNS2=%f, df2_DNS2=%f, f1_DNS2=%f, f2_DNS2=%f, weight_DNS2=%f\n",df1_NLO_DNS2,df2_NLO_DNS2,f1_NLO,f2_NLO,weight_NLO_DNS2);
565  }
566 
567 
568  fillPythiaEvent(mEvent);
569 
570  }
571 
572  return kStOK;
573 }
574 
575 void StMCAsymMaker::fillPythiaEvent(StPythiaEvent* pythia)
576 {
577  g2t_event_st* eventTable = Pg2t_event->GetTable();
578 
579  pythia->setRunId(eventTable->n_run);
580  pythia->setEventId(eventTable->n_event);
581 
582  if (mcEvent && mcEvent->primaryVertex())
583  pythia->setVertex(mcEvent->primaryVertex()->position().xyz());
584 
585  g2t_pythia_st* pythiaTable = Pg2t_pythia->GetTable();
586 
587  pythia->setProcessId(pythiaTable->subprocess_id);
588  pythia->setS(pythiaTable->mand_s);
589  pythia->setT(pythiaTable->mand_t);
590  pythia->setU(pythiaTable->mand_u);
591  pythia->setPt(pythiaTable->hard_p);
592  pythia->setCosTheta(pythiaTable->cos_th);
593  pythia->setX1(pythiaTable->bjor_1);
594  pythia->setX2(pythiaTable->bjor_2);
595  pythia->setMstu72(pythiaTable->mstu72);
596  pythia->setMstu73(pythiaTable->mstu73);
597  pythia->setMstp111(pythiaTable->mstp111);
598  pythia->setPartonALL(partonic_all);
599 
600  pythia->setDF1(StPythiaEvent::LO, df1_LO);
601  pythia->setDF1(StPythiaEvent::NLO, df1_NLO);
602  pythia->setDF1(StPythiaEvent::ZERO, df1_NLO_g0);
603  pythia->setDF1(StPythiaEvent::MAX, df1_NLO_gmax);
604  pythia->setDF1(StPythiaEvent::MIN, df1_NLO_gmin);
605  pythia->setDF1(StPythiaEvent::M015, df1_NLO_m015);
606  pythia->setDF1(StPythiaEvent::M030, df1_NLO_m030);
607  pythia->setDF1(StPythiaEvent::M045, df1_NLO_m045);
608  pythia->setDF1(StPythiaEvent::M060, df1_NLO_m060);
609  pythia->setDF1(StPythiaEvent::M075, df1_NLO_m075);
610  pythia->setDF1(StPythiaEvent::M090, df1_NLO_m090);
611  pythia->setDF1(StPythiaEvent::M105, df1_NLO_m105);
612  pythia->setDF1(StPythiaEvent::P030, df1_NLO_p030);
613  pythia->setDF1(StPythiaEvent::P045, df1_NLO_p045);
614  pythia->setDF1(StPythiaEvent::P060, df1_NLO_p060);
615  pythia->setDF1(StPythiaEvent::P070, df1_NLO_p070);
616  pythia->setDF1(StPythiaEvent::GS_NLOA, df1_NLO_GSA);
617  pythia->setDF1(StPythiaEvent::GS_NLOB, df1_NLO_GSB);
618  pythia->setDF1(StPythiaEvent::GS_NLOC, df1_NLO_GSC);
619  pythia->setDF1(StPythiaEvent::DSSV, df1_NLO_DSSV);
620  pythia->setDF1(StPythiaEvent::DSSV2009a, df1_NLO_DSSV2009a);
621  pythia->setDF1(StPythiaEvent::LSS1, df1_NLO_LSS1);
622  pythia->setDF1(StPythiaEvent::LSS2, df1_NLO_LSS2);
623  pythia->setDF1(StPythiaEvent::LSS3, df1_NLO_LSS3);
624  pythia->setDF1(StPythiaEvent::LSS2010_delGpos, df1_NLO_LSS2010_delGpos);
625  pythia->setDF1(StPythiaEvent::LSS2010_chsign_delG, df1_NLO_LSS2010_chsign_delG);
626  pythia->setDF1(StPythiaEvent::AAC1, df1_NLO_AAC1);
627  pythia->setDF1(StPythiaEvent::AAC2, df1_NLO_AAC2);
628  pythia->setDF1(StPythiaEvent::AAC3, df1_NLO_AAC3);
629  pythia->setDF1(StPythiaEvent::BB1, df1_NLO_BB1);
630  pythia->setDF1(StPythiaEvent::BB2, df1_NLO_BB2);
631  pythia->setDF1(StPythiaEvent::BB2010, df1_NLO_BB2010);
632  pythia->setDF1(StPythiaEvent::DNS1, df1_NLO_DNS1);
633  pythia->setDF1(StPythiaEvent::DNS2, df1_NLO_DNS2);
634 
635  pythia->setDF2(StPythiaEvent::LO, df2_LO);
636  pythia->setDF2(StPythiaEvent::NLO, df2_NLO);
637  pythia->setDF2(StPythiaEvent::ZERO, df2_NLO_g0);
638  pythia->setDF2(StPythiaEvent::MAX, df2_NLO_gmax);
639  pythia->setDF2(StPythiaEvent::MIN, df2_NLO_gmin);
640  pythia->setDF2(StPythiaEvent::M015, df2_NLO_m015);
641  pythia->setDF2(StPythiaEvent::M030, df2_NLO_m030);
642  pythia->setDF2(StPythiaEvent::M045, df2_NLO_m045);
643  pythia->setDF2(StPythiaEvent::M060, df2_NLO_m060);
644  pythia->setDF2(StPythiaEvent::M075, df2_NLO_m075);
645  pythia->setDF2(StPythiaEvent::M090, df2_NLO_m090);
646  pythia->setDF2(StPythiaEvent::M105, df2_NLO_m105);
647  pythia->setDF2(StPythiaEvent::P030, df2_NLO_p030);
648  pythia->setDF2(StPythiaEvent::P045, df2_NLO_p045);
649  pythia->setDF2(StPythiaEvent::P060, df2_NLO_p060);
650  pythia->setDF2(StPythiaEvent::P070, df2_NLO_p070);
651  pythia->setDF2(StPythiaEvent::GS_NLOA, df2_NLO_GSA);
652  pythia->setDF2(StPythiaEvent::GS_NLOB, df2_NLO_GSB);
653  pythia->setDF2(StPythiaEvent::GS_NLOC, df2_NLO_GSC);
654  pythia->setDF2(StPythiaEvent::DSSV, df2_NLO_DSSV);
655  pythia->setDF2(StPythiaEvent::DSSV2009a, df2_NLO_DSSV2009a);
656  pythia->setDF2(StPythiaEvent::LSS1, df2_NLO_LSS1);
657  pythia->setDF2(StPythiaEvent::LSS2, df2_NLO_LSS2);
658  pythia->setDF2(StPythiaEvent::LSS3, df2_NLO_LSS3);
659  pythia->setDF2(StPythiaEvent::LSS2010_delGpos, df2_NLO_LSS2010_delGpos);
660  pythia->setDF2(StPythiaEvent::LSS2010_chsign_delG, df2_NLO_LSS2010_chsign_delG);
661  pythia->setDF2(StPythiaEvent::AAC1, df2_NLO_AAC1);
662  pythia->setDF2(StPythiaEvent::AAC2, df2_NLO_AAC2);
663  pythia->setDF2(StPythiaEvent::AAC3, df2_NLO_AAC3);
664  pythia->setDF2(StPythiaEvent::BB1, df2_NLO_BB1);
665  pythia->setDF2(StPythiaEvent::BB2, df2_NLO_BB2);
666  pythia->setDF2(StPythiaEvent::BB2010, df2_NLO_BB2010);
667  pythia->setDF2(StPythiaEvent::DNS1, df2_NLO_DNS1);
668  pythia->setDF2(StPythiaEvent::DNS2, df2_NLO_DNS2);
669 
670  pythia->setF1(StPythiaEvent::LO, f1_LO);
671  pythia->setF1(StPythiaEvent::NLO, f1_NLO);
672 
673  pythia->setF2(StPythiaEvent::LO, f2_LO);
674  pythia->setF2(StPythiaEvent::NLO, f2_NLO);
675 
676  particle_st* particleTable = particleTabPtr->GetTable();
677 
678  for (int i = 0; i < particleTabPtr->GetNRows(); ++i)
679  pythia->addParticle(TParticle(particleTable[i].idhep,
680  particleTable[i].isthep,
681  particleTable[i].jmohep[0],
682  particleTable[i].jmohep[1],
683  particleTable[i].jdahep[0],
684  particleTable[i].jdahep[1],
685  TLorentzVector(particleTable[i].phep),
686  TLorentzVector(particleTable[i].vhep)));
687 
688 }
689 
690 //Gehrmann-Stirling NL0 Set A
691 Double_t StMCAsymMaker::get_polPDF_NLO_GSA(int flavor, double x, double Q2){
692 
693  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
694  double pdf=1000;
695  int polset_NLO_GSA=201;
696  int polid=0;
697 
698  // cout<<"get_polPDF_NLO_GSA: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
699  if ((Q2>=1.0)&&(Q2<=1E+06)) polar_(&polset_NLO_GSA, &x, &Q2, parpol, &polid);
700  // cout<<"get_polPDF_NLO_GSA: UV="<<parpol[0]<<" DV="<<parpol[1]<<" UB="<<parpol[3]<<" DB="<<parpol[4]<<" S+Sbar="<<parpol[5]<<" GL="<<parpol[2]<<endl;
701 
702 
703  //parpol[0] = uv - ubar so it is only the valence contribution. Need to add back in ubar for pythia since we don't know the difference between the sea and valence
704  if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
705  //parpol[1] = dv - dbar
706  if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
707  //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
708  if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark
709  //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
710  if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
711  //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
712  if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
713  if (flavor==21) pdf=parpol[5]; //gluon
714  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
715 
716  return pdf;
717 }
718 
719 //Gehrmann-Stirling NL0 Set B
720 Double_t StMCAsymMaker::get_polPDF_NLO_GSB(int flavor, double x, double Q2){
721 
722  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
723  double pdf=1000;
724  int polset_NLO_GSB=202;
725  int polid=0;
726 
727  //cout<<"get_polPDF_NLO_GSB: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
728  if ((Q2>=1.0)&&(Q2<=1E+06)) polar_(&polset_NLO_GSB, &x, &Q2, parpol, &polid);
729  //cout<<"get_polPDF_NLO_GSB: UV="<<parpol[0]<<" DV="<<parpol[1]<<" UB="<<parpol[3]<<" DB="<<parpol[4]<<" S+Sbar="<<parpol[5]<<" GL="<<parpol[2]<<endl;
730 
731 
732  //parpol[0] = uv - ubar so it is only the valence contribution. Need to add back in ubar for pythia since we don't know the difference between the sea and valence
733  if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
734  //parpol[1] = dv - dbar
735  if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
736  //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
737  if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark
738  //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
739  if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
740  //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
741  if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
742  if (flavor==21) pdf=parpol[5]; //gluon
743  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
744 
745  return pdf;
746 }
747 
748 //Gehrmann-Stirling NL0 Set C
749 Double_t StMCAsymMaker::get_polPDF_NLO_GSC(int flavor, double x, double Q2){
750 
751  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
752  double pdf=1000;
753  int polset_NLO_GSC=203;
754  int polid=0;
755 
756  //cout<<"get_polPDF_NLO_GSC: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
757  if ((Q2>=1.0)&&(Q2<=1E+06)) polar_(&polset_NLO_GSC, &x, &Q2, parpol, &polid);
758  //cout<<"get_polPDF_NLO_GSC: UV="<<parpol[0]<<" DV="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" S+Sbar="<<parpol[4]<<" GL="<<parpol[5]<<endl;
759 
760  //parpol[0] = uv - ubar so it is only the valence contribution. Need to add back in ubar for pythia since we don't know the difference between the sea and valence
761  if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
762  //parpol[1] = dv - dbar
763  if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
764  //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
765  if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark
766  //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
767  if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
768  //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
769  if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
770  if (flavor==21) pdf=parpol[5]; //gluon
771  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
772 
773 
774  return pdf;
775 }
776 
777 
778 
779 
780 
781 //GRSV LO standard
782 Double_t StMCAsymMaker::get_polPDF_LO(int flavor, double x, double Q2){
783 
784  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
785  double pdf=1000;
786  int polset_LO=101;
787  int polid=0;
788  //cout<<"get_polPDF_LO: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
789 
790  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_LO, &x, &Q2, parpol, &polid);
791  //cout <<"getpolPDF_LO: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
792 
793  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
794  if (flavor==2) pdf=parpol[0]; //uv + usea quark
795  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
796  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
797  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
798  if (flavor==21) pdf=parpol[5]; //gluon
799  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
800 
801  return pdf;
802 }
803 
804 
805 //GRSV NLO standard
806 Double_t StMCAsymMaker::get_polPDF_NLO(int flavor, double x, double Q2){
807 
808  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
809  double pdf=1000;
810  int polset_NLO=102;
811  int polid=0;
812  //cout<<"get_polPDF_NLO: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
813 
814  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO, &x, &Q2, parpol, &polid);
815  //cout <<"get_polPDF_NLO: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
816 
817  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
818  if (flavor==2) pdf=parpol[0]; //uv + usea quark
819  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
820  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
821  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
822  if (flavor==21) pdf=parpol[5]; //gluon
823  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
824 
825  return pdf;
826 }
827 
828 //GRSV NLO g0
829 Double_t StMCAsymMaker::get_polPDF_NLO_g0(int flavor, double x, double Q2){
830 
831  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
832  double pdf=1000;
833  int polset_NLO_g0=103;
834  int polid=0;
835  //cout<<"get_polPDF_NLO_g0: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
836 
837  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_g0, &x, &Q2, parpol, &polid);
838  //cout <<"getpolPDF_NLO_g0: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
839 
840  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
841  if (flavor==2) pdf=parpol[0]; //uv + usea quark
842  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
843  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
844  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
845  if (flavor==21) pdf=parpol[5]; //gluon
846  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
847 
848  return pdf;
849 }
850 
851 
852 //GRSV NLO gmax
853 Double_t StMCAsymMaker::get_polPDF_NLO_gmax(int flavor, double x, double Q2){
854 
855  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
856  double pdf=1000;
857  int polset_NLO_gmax=104;
858  int polid=0;
859  //cout<<"get_polPDF_NLO_gmax: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
860 
861  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_gmax, &x, &Q2, parpol, &polid);
862  //cout <<"getpolPDF_NLO_gmax: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
863 
864  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
865  if (flavor==2) pdf=parpol[0]; //uv + usea quark
866  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
867  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
868  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
869  if (flavor==21) pdf=parpol[5]; //gluon
870  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
871 
872  return pdf;
873 }
874 
875 
876 //GRSV NLO gmin
877 Double_t StMCAsymMaker::get_polPDF_NLO_gmin(int flavor, double x, double Q2){
878 
879  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
880  double pdf=1000;
881  int polset_NLO_gmin=105;
882  int polid=0;
883  //cout<<"get_polPDF_NLO_gmin: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
884 
885  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_gmin, &x, &Q2, parpol, &polid);
886  //cout <<"getpolPDF_NLO_gmin: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
887 
888  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
889  if (flavor==2) pdf=parpol[0]; //uv + usea quark
890  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
891  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
892  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
893  if (flavor==21) pdf=parpol[5]; //gluon
894  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
895 
896  return pdf;
897 }
898 
899 //GRSV NLO m015
900 Double_t StMCAsymMaker::get_polPDF_NLO_m015(int flavor, double x, double Q2){
901 
902  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
903  double pdf=1000;
904  int polset_NLO_m015=106;
905  int polid=0;
906  //cout<<"get_polPDF_NLO_m015: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
907 
908  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m015, &x, &Q2, parpol, &polid);
909  //cout <<"getpolPDF_NLO_m015: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
910 
911  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
912  if (flavor==2) pdf=parpol[0]; //uv + usea quark
913  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
914  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
915  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
916  if (flavor==21) pdf=parpol[5]; //gluon
917  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
918 
919  return pdf;
920 }
921 
922 //GRSV NLO m030
923 Double_t StMCAsymMaker::get_polPDF_NLO_m030(int flavor, double x, double Q2){
924 
925  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
926  double pdf=1000;
927  int polset_NLO_m030=107;
928  int polid=0;
929  //cout<<"get_polPDF_NLO_m030: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
930 
931  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m030, &x, &Q2, parpol, &polid);
932  //cout <<"getpolPDF_NLO_m030: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
933 
934  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
935  if (flavor==2) pdf=parpol[0]; //uv + usea quark
936  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
937  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
938  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
939  if (flavor==21) pdf=parpol[5]; //gluon
940  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
941 
942  return pdf;
943 }
944 
945 //GRSV NLO m045
946 Double_t StMCAsymMaker::get_polPDF_NLO_m045(int flavor, double x, double Q2){
947 
948  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
949  double pdf=1000;
950  int polset_NLO_m045=108;
951  int polid=0;
952  //cout<<"get_polPDF_NLO_m045: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
953 
954  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m045, &x, &Q2, parpol, &polid);
955  //cout <<"getpolPDF_NLO_m045: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
956 
957  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
958  if (flavor==2) pdf=parpol[0]; //uv + usea quark
959  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
960  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
961  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
962  if (flavor==21) pdf=parpol[5]; //gluon
963  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
964 
965  return pdf;
966 }
967 
968 //GRSV NLO m060
969 Double_t StMCAsymMaker::get_polPDF_NLO_m060(int flavor, double x, double Q2){
970 
971  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
972  double pdf=1000;
973  int polset_NLO_m060=109;
974  int polid=0;
975  //cout<<"get_polPDF_NLO_m060: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
976 
977  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m060, &x, &Q2, parpol, &polid);
978  //cout <<"getpolPDF_NLO_m060: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
979 
980  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
981  if (flavor==2) pdf=parpol[0]; //uv + usea quark
982  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
983  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
984  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
985  if (flavor==21) pdf=parpol[5]; //gluon
986  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
987 
988  return pdf;
989 }
990 
991 //GRSV NLO m075
992 Double_t StMCAsymMaker::get_polPDF_NLO_m075(int flavor, double x, double Q2){
993 
994  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
995  double pdf=1000;
996  int polset_NLO_m075=110;
997  int polid=0;
998  //cout<<"get_polPDF_NLO_m075: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
999 
1000  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m075, &x, &Q2, parpol, &polid);
1001  //cout <<"getpolPDF_NLO_m075: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1002 
1003  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1004  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1005  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1006  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1007  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1008  if (flavor==21) pdf=parpol[5]; //gluon
1009  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1010 
1011  return pdf;
1012 }
1013 
1014 //GRSV NLO m090
1015 Double_t StMCAsymMaker::get_polPDF_NLO_m090(int flavor, double x, double Q2){
1016 
1017  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1018  double pdf=1000;
1019  int polset_NLO_m090=111;
1020  int polid=0;
1021  //cout<<"get_polPDF_NLO_m090: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1022 
1023  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m090, &x, &Q2, parpol, &polid);
1024  //cout <<"getpolPDF_NLO_m090: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1025 
1026  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1027  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1028  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1029  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1030  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1031  if (flavor==21) pdf=parpol[5]; //gluon
1032  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1033 
1034  return pdf;
1035 }
1036 
1037 //GRSV NLO m105
1038 Double_t StMCAsymMaker::get_polPDF_NLO_m105(int flavor, double x, double Q2){
1039 
1040  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1041  double pdf=1000;
1042  int polset_NLO_m105=112;
1043  int polid=0;
1044  //cout<<"get_polPDF_NLO_m105: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1045 
1046  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m105, &x, &Q2, parpol, &polid);
1047  //cout <<"getpolPDF_NLO_m105: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1048 
1049  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1050  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1051  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1052  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1053  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1054  if (flavor==21) pdf=parpol[5]; //gluon
1055  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1056 
1057  return pdf;
1058 }
1059 
1060 //GRSV NLO p030
1061 Double_t StMCAsymMaker::get_polPDF_NLO_p030(int flavor, double x, double Q2){
1062 
1063  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1064  double pdf=1000;
1065  int polset_NLO_p030=113;
1066  int polid=0;
1067  //cout<<"get_polPDF_NLO_p030: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1068 
1069  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p030, &x, &Q2, parpol, &polid);
1070  //cout <<"getpolPDF_NLO_p030: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1071 
1072  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1073  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1074  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1075  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1076  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1077  if (flavor==21) pdf=parpol[5]; //gluon
1078  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1079 
1080  return pdf;
1081 }
1082 
1083 //GRSV NLO p045
1084 Double_t StMCAsymMaker::get_polPDF_NLO_p045(int flavor, double x, double Q2){
1085 
1086  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1087  double pdf=1000;
1088  int polset_NLO_p045=114;
1089  int polid=0;
1090  //cout<<"get_polPDF_NLO_p045: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1091 
1092  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p045, &x, &Q2, parpol, &polid);
1093  //cout <<"getpolPDF_NLO_p045: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1094 
1095  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1096  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1097  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1098  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1099  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1100  if (flavor==21) pdf=parpol[5]; //gluon
1101  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1102 
1103  return pdf;
1104 }
1105 
1106 //GRSV NLO p060
1107 Double_t StMCAsymMaker::get_polPDF_NLO_p060(int flavor, double x, double Q2){
1108 
1109  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1110  double pdf=1000;
1111  int polset_NLO_p060=115;
1112  int polid=0;
1113  //cout<<"get_polPDF_NLO_p060: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1114 
1115  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p060, &x, &Q2, parpol, &polid);
1116  //cout <<"getpolPDF_NLO_p060: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1117 
1118  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1119  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1120  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1121  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1122  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1123  if (flavor==21) pdf=parpol[5]; //gluon
1124  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1125 
1126  return pdf;
1127 }
1128 
1129 //GRSV NLO p070
1130 Double_t StMCAsymMaker::get_polPDF_NLO_p070(int flavor, double x, double Q2){
1131 
1132  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1133  double pdf=1000;
1134  int polset_NLO_p070=116;
1135  int polid=0;
1136  //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1137 
1138  if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p070, &x, &Q2, parpol, &polid);
1139  //cout <<"getpolPDF_NLO_p070: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1140 
1141  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1142  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1143  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1144  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1145  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1146  if (flavor==21) pdf=parpol[5]; //gluon
1147  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1148 
1149  return pdf;
1150 }
1151 
1152 //DSSV NLO
1153 Double_t StMCAsymMaker::get_polPDF_NLO_DSSV(int flavor, double x, double Q2){
1154 
1155  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1156  double pdf=1000;
1157  int polset_NLO_DSSV=301;
1158  int polid=0;
1159  //cout<<"get_polPDF_NLO_DSSV: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1160 
1161  if (Q2>=1.0&&Q2<=1.0e5) polar_(&polset_NLO_DSSV, &x, &Q2, parpol, &polid);
1162  //cout <<"getpolPDF_NLO_DSSV: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1163 
1164  if (flavor==1) pdf=parpol[1]+parpol[2]; //dv + dsea quark
1165  if (flavor==2) pdf=parpol[0]+parpol[3]; //uv + usea quark
1166  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1167  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1168  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1169  if (flavor==21) pdf=parpol[5]; //gluon
1170  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1171 
1172  return pdf;
1173 }
1174 
1175 //DSSV 2009a NLO
1176 Double_t StMCAsymMaker::get_polPDF_NLO_DSSV2009a(int flavor, double x, double Q2){
1177 
1178  double duv=0,ddv=0,dubar=0,ddbar=0,dstr=0,dglu=0;
1179 
1180  if (x>=1.0e-5&&x<=1.0&&Q2>=1.0&&Q2<=1.0e5) {
1181  dssvfit2009a_(&x,&Q2,&duv,&ddv,&dubar,&ddbar,&dstr,&dglu);
1182  switch (flavor) {
1183  case 1: return (ddv+ddbar)/x; //dv + dsea quark
1184  case 2: return (duv+dubar)/x; //uv + usea quark
1185  case -1: return ddbar/x; //dbar==dsea quark
1186  case -2: return dubar/x; //ubar==usea quark
1187  case 3:
1188  case -3: return dstr/x; //s==sbar quark
1189  case 21: return dglu/x; //gluon
1190  case 4:
1191  case -4:
1192  case 5:
1193  case -5:
1194  case 6:
1195  case -6: return dstr/x;
1196  }
1197  }
1198  return 1000;
1199 }
1200 
1201 //LSS NLO Scenario 1
1202 Double_t StMCAsymMaker::get_polPDF_NLO_LSS1(int flavor, double x, double Q2){
1203 
1204  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1205  double pdf=1000;
1206  int polset_NLO_LSS1=401;
1207  int polid=0;
1208  // cout<<"get_polPDF_NLO_LSS1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1209 
1210  if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS1, &x, &Q2, parpol, &polid);
1211  // cout <<"getpolPDF_NLO_LSS1: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1212 
1213  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1214  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1215  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1216  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1217  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1218  if (flavor==21) pdf=parpol[5]; //gluon
1219  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1220 
1221  return pdf;
1222 }
1223 
1224 //LSS NLO Scenario 2
1225 Double_t StMCAsymMaker::get_polPDF_NLO_LSS2(int flavor, double x, double Q2){
1226 
1227  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1228  double pdf=1000;
1229  int polset_NLO_LSS2=402;
1230  int polid=0;
1231  // cout<<"get_polPDF_NLO_LSS2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1232 
1233  if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS2, &x, &Q2, parpol, &polid);
1234  // cout <<"getpolPDF_NLO_LSS2: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1235 
1236  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1237  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1238  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1239  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1240  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1241  if (flavor==21) pdf=parpol[5]; //gluon
1242  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1243 
1244  return pdf;
1245 }
1246 
1247 //LSS NLO Scenario 3
1248 Double_t StMCAsymMaker::get_polPDF_NLO_LSS3(int flavor, double x, double Q2){
1249 
1250  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1251  double pdf=1000;
1252  int polset_NLO_LSS3=403;
1253  int polid=0;
1254  // cout<<"get_polPDF_NLO_LSS3: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1255 
1256  if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS3, &x, &Q2, parpol, &polid);
1257  // cout <<"getpolPDF_NLO_LSS3: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1258 
1259  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1260  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1261  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1262  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1263  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1264  if (flavor==21) pdf=parpol[5]; //gluon
1265  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1266 
1267  return pdf;
1268 }
1269 
1270 //LSS2010 pos x*DeltaG
1271 Double_t StMCAsymMaker::get_polPDF_NLO_LSS2010_delGpos(int flavor, double x, double Q2)
1272 {
1273  static int iini = 0;
1274  int iset = 1;
1275  double uub, ddb, u, d, ub, db, st, gl;
1276 
1277  if (iini == 0) intini_.iini = 0;
1278 
1279  lss2010_(&iset,&x,&Q2,&uub,&ddb,&u,&d,&ub,&db,&st,&gl);
1280 
1281  if (iini == 0) iini = 1;
1282 
1283  switch (flavor) {
1284  case 1: return ddb/x;
1285  case 2: return uub/x;
1286  case -1: return db/x;
1287  case -2: return ub/x;
1288  case 3:
1289  case 4:
1290  case 5:
1291  case 6:
1292  case -3:
1293  case -4:
1294  case -5:
1295  case -6: return st/x;
1296  case 21: return gl/x;
1297  }
1298 
1299  return 1000;
1300 }
1301 
1302 //LSS2010 node x*DeltaG
1303 Double_t StMCAsymMaker::get_polPDF_NLO_LSS2010_chsign_delG(int flavor, double x, double Q2)
1304 {
1305  static int iini = 0;
1306  int iset = 2;
1307  double uub, ddb, u, d, ub, db, st, gl;
1308 
1309  if (iini == 0) intini_.iini = 0;
1310 
1311  lss2010_(&iset,&x,&Q2,&uub,&ddb,&u,&d,&ub,&db,&st,&gl);
1312 
1313  if (iini == 0) iini = 1;
1314 
1315  switch (flavor) {
1316  case 1: return ddb/x;
1317  case 2: return uub/x;
1318  case -1: return db/x;
1319  case -2: return ub/x;
1320  case 3:
1321  case 4:
1322  case 5:
1323  case 6:
1324  case -3:
1325  case -4:
1326  case -5:
1327  case -6: return st/x;
1328  case 21: return gl/x;
1329  }
1330 
1331  return 1000;
1332 }
1333 
1334 //AAC NLO Scenario 1
1335 Double_t StMCAsymMaker::get_polPDF_NLO_AAC1(int flavor, double x, double Q2){
1336 
1337  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1338  double pdf=1000;
1339  int polset_NLO_AAC1=501;
1340  int polid=0;
1341  //cout<<"get_polPDF_NLO_AAC1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1342 
1343  if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC1, &x, &Q2, parpol, &polid);
1344  //cout <<"getpolPDF_NLO_AAC1: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1345 
1346  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1347  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1348  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1349  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1350  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1351  if (flavor==21) pdf=parpol[5]; //gluon
1352  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1353 
1354  return pdf;
1355 }
1356 
1357 //AAC NLO Scenario 2
1358 Double_t StMCAsymMaker::get_polPDF_NLO_AAC2(int flavor, double x, double Q2){
1359 
1360  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1361  double pdf=1000;
1362  int polset_NLO_AAC2=502;
1363  int polid=0;
1364  //cout<<"get_polPDF_NLO_AAC2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1365 
1366  if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC2, &x, &Q2, parpol, &polid);
1367  //cout <<"getpolPDF_NLO_AAC2: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1368 
1369  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1370  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1371  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1372  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1373  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1374  if (flavor==21) pdf=parpol[5]; //gluon
1375  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1376 
1377  return pdf;
1378 }
1379 
1380 //AAC NLO Scenario 3
1381 Double_t StMCAsymMaker::get_polPDF_NLO_AAC3(int flavor, double x, double Q2){
1382 
1383  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1384  double pdf=1000;
1385  int polset_NLO_AAC3=503;
1386  int polid=0;
1387  //cout<<"get_polPDF_NLO_AAC3: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1388 
1389  if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC3, &x, &Q2, parpol, &polid);
1390  //cout <<"getpolPDF_NLO_AAC3: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1391 
1392  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1393  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1394  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1395  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1396  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1397  if (flavor==21) pdf=parpol[5]; //gluon
1398  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1399 
1400  return pdf;
1401 }
1402 
1403 //BB NLO Scenario 1
1404 Double_t StMCAsymMaker::get_polPDF_NLO_BB1(int flavor, double x, double Q2){
1405 
1406  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1407  double pdf=1000;
1408  int polset_NLO_BB1=601;
1409  int polid=0;
1410  //cout<<"get_polPDF_NLO_BB1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1411 
1412  if (Q2>=1.0&&Q2<=1.0e6) polar_(&polset_NLO_BB1, &x, &Q2, parpol, &polid);
1413  //cout <<"get_polPDF_NLO_BB1: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1414 
1415  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1416  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1417  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1418  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1419  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1420  if (flavor==21) pdf=parpol[5]; //gluon
1421  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1422 
1423  return pdf;
1424 }
1425 
1426 //BB NLO Scenario 2
1427 Double_t StMCAsymMaker::get_polPDF_NLO_BB2(int flavor, double x, double Q2){
1428 
1429  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1430  double pdf=1000;
1431  int polset_NLO_BB2=602;
1432  int polid=0;
1433  //cout<<"get_polPDF_NLO_BB2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1434 
1435  if (Q2>=1.0&&Q2<=1.0e6) polar_(&polset_NLO_BB2, &x, &Q2, parpol, &polid);
1436  //cout <<"get_polPDF_NLO_BB2: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1437 
1438  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1439  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1440  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1441  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1442  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1443  if (flavor==21) pdf=parpol[5]; //gluon
1444  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1445 
1446  return pdf;
1447 }
1448 
1449 //BB NL0 2010
1450 Double_t StMCAsymMaker::get_polPDF_NLO_BB2010(int flavor, double x, double Q2)
1451 {
1452  static int iini = 0;
1453  int iset = 1;
1454  double uv, duv, dv, ddv, gl, dgl, qb, dqb, g1p, dg1p, g1n, dg1n;
1455 
1456  if (iini == 0) intini_.iini = 0;
1457 
1458  polpdf_(&iset,&x,&Q2,&uv,&duv,&dv,&ddv,&gl,&dgl,&qb,&dqb,&g1p,&dg1p,&g1n,&dg1n);
1459 
1460  if (iini == 0) iini = 1;
1461 
1462  switch (flavor) {
1463  case 1: return (dv+qb)/x;
1464  case 2: return (uv+qb)/x;
1465  case 3:
1466  case 4:
1467  case 5:
1468  case 6:
1469  case -1:
1470  case -2:
1471  case -3:
1472  case -4:
1473  case -5:
1474  case -6: return qb/x;
1475  case 21: return gl/x;
1476  }
1477 
1478  return 1000;
1479 }
1480 
1481 //DNS NLO Scenario 1
1482 Double_t StMCAsymMaker::get_polPDF_NLO_DNS1(int flavor, double x, double Q2){
1483 
1484  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1485  double pdf=1000;
1486  int polset_NLO_DNS1=701;
1487  int polid=0;
1488  //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1489 
1490  if (Q2>=1.0&&Q2<=5.0e4) polar_(&polset_NLO_DNS1, &x, &Q2, parpol, &polid);
1491  //cout <<"getpolPDF_NLO_p070: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1492 
1493  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1494  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1495  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1496  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1497  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1498  if (flavor==21) pdf=parpol[5]; //gluon
1499  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1500 
1501  return pdf;
1502 }
1503 
1504 //DNS NLO Scenario 2
1505 Double_t StMCAsymMaker::get_polPDF_NLO_DNS2(int flavor, double x, double Q2){
1506 
1507  double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
1508  double pdf=1000;
1509  int polset_NLO_DNS2=702;
1510  int polid=0;
1511  //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
1512 
1513  if (Q2>=1.0&&Q2<=5.0e4) polar_(&polset_NLO_DNS2, &x, &Q2, parpol, &polid);
1514  //cout <<"getpolPDF_NLO_p070: U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
1515 
1516  if (flavor==1) pdf=parpol[1]; //dv + dsea quark
1517  if (flavor==2) pdf=parpol[0]; //uv + usea quark
1518  if (flavor==-1) pdf=parpol[3]; //dbar==dsea quark
1519  if (flavor==-2) pdf=parpol[2]; //ubar==usea quark
1520  if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
1521  if (flavor==21) pdf=parpol[5]; //gluon
1522  if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
1523 
1524  return pdf;
1525 }
1526 
1527 //Returns unpolarized CTEQ
1528 Double_t StMCAsymMaker::get_unpolPDF_LO(int flavor, double x, double Q2){
1529 
1530  Double_t pdf=0.0;
1531  Int_t iset=3;//LO
1532  Int_t er=0;
1533  Int_t fl=10;
1534 
1535  if (flavor==1) fl=2;
1536  if (flavor==2) fl=1;
1537  if (flavor==-1) fl=-2;
1538  if (flavor==-2) fl=-1;
1539  if (flavor==21) fl=0;
1540  if (flavor==3) fl=3;
1541  if (flavor==-3) fl=-3;
1542  if (flavor==4) fl=4;
1543  if (flavor==-4) fl=-4;
1544  if (flavor==5) fl=5;
1545  if (flavor==-5) fl=-5;
1546 
1547  double Q=pow(Q2,0.5);
1548  pdf=ctq5pd_(&iset,&fl,&x,&Q,&er);
1549  if (er!=0) pdf=0.0;
1550 
1551  return pdf;
1552 }
1553 
1554 Double_t StMCAsymMaker::get_unpolPDF_NLO(int flavor, double x, double Q2){
1555 
1556  Double_t pdf=0.0;
1557  Int_t iset=1;//NLO MSbar scheme
1558  Int_t er=0;
1559  Int_t fl=10;
1560 
1561  if (flavor==1) fl=2;
1562  if (flavor==2) fl=1;
1563  if (flavor==-1) fl=-2;
1564  if (flavor==-2) fl=-1;
1565  if (flavor==21) fl=0;
1566  if (flavor==3) fl=3;
1567  if (flavor==-3) fl=-3;
1568  if (flavor==4) fl=4;
1569  if (flavor==-4) fl=-4;
1570  if (flavor==5) fl=5;
1571  if (flavor==-5) fl=-5;
1572 
1573  double Q=pow(Q2,0.5);
1574  pdf=ctq5pd_(&iset,&fl,&x,&Q,&er);
1575  if (er!=0) pdf=0.0;
1576 
1577  return pdf;
1578 }
1579 
1580 
1581 Double_t StMCAsymMaker::getPartonicALL(double s, double t, double u, int sub, int inA, int inB, int outA, int outB){
1582 
1583  //Werner definitions:
1584  //1: qq'->qq' (qqbar'->qqbar') 2: qq->qq 3: qqbar->q'qbar' 4: qqbar->qqbar 5: qqbar->gg 6: gg->qqbar 7: qg->qg 8: gg->gg
1585  //PYTHIA definitions:
1586  //1: 11a 2:11b 3: 12a 4:11 and 12b 5: 13 6: 53 7: 28 8: 68
1587  //NOTES:
1588  // 3==5==6==-1 1==7 1!=2 and 1!=4
1589 
1590  double N1,N2,N3,N4,N5,N6,N7,N8;
1591  double D1,D2,D3,D4,D5,D6,D7,D8;
1592  double all=-10;
1593 
1594  num_(&s,&t,&u,&N1,&N2,&N3,&N4,&N5,&N6,&N7,&N8);
1595  denom_(&s,&t,&u,&D1,&D2,&D3,&D4,&D5,&D6,&D7,&D8);
1596  if (0){
1597  cout<<"s="<<s<<" t="<<t<<" u="<<u<<" sub="<<sub<<" inA="<<inA<<" inB="<<inB<<" outA="<<outA<<" outB="<<outB<<endl;
1598  cout<<" 1="<<N1<<" "<<D1<<endl;
1599  cout<<" 2="<<N2<<" "<<D2<<endl;
1600  cout<<" 3="<<N3<<" "<<D3<<endl;
1601  cout<<" 4="<<N4<<" "<<D4<<endl;
1602  cout<<" 5="<<N5<<" "<<D5<<endl;
1603  cout<<" 6="<<N6<<" "<<D6<<endl;
1604  cout<<" 7="<<N7<<" "<<D7<<endl;
1605  cout<<" 8="<<N8<<" "<<D8<<endl;
1606  }
1607 
1608 
1609  if ((sub==11)&&(abs(inA)!=abs(inB))) all=N1/D1;
1610  if ((sub==11)&&(inA==inB)&&(outA==outB)) all=N2/D2;
1611  //This line added as bug fix 07/17/08 RHF
1612  //Before then these events given all=N2/D2
1613  if ((sub==11)&&(inA==(-1*inB))&&(outA==(-1*outB))&&(inA==outA)&&(inB==outB)) all=N4/D4;
1614  if ((sub==12)&&(abs(inA)!=abs(outA))) all=N3/D3;
1615  if ((sub==12)&&(abs(inA)==abs(outA))) all=N4/D4;
1616  if (sub==13) all=N5/D5;
1617  if (sub==53) all=N6/D6;
1618  if (sub==28) all=N7/D7;
1619  if (sub==68) all=N8/D8;
1620 
1621  //prompt photon subprocesses
1622  if (sub==29) all=N7/D7;//q_i + g -> q_i + gamma
1623  if (sub==14) all=N6/D6;//q_i+qbar_i -> g + gamma
1624  if (sub==18) all=N6/D6;//q_i+qbar_i -> gamma + gamma
1625 
1626  return all;
1627 }
1628 
1629 Double_t StMCAsymMaker::getProtonA1(Double_t x,Double_t Q2){
1630 
1631  // This piece coded by Pibero
1632  // Virtual photon-nucleon asymmetry measured from deep inelatic scattering
1633  // of polarized charged lepton beams from polarized targets. It is defined
1634  // as the sum of polarized PDF weighted by charge square over the sum of
1635  // unpolarized PDF weighted by charge square for all quark/antiquark flavors.
1636 
1637  // Flavors are PDG codes for the quarks and gluons:
1638  // down=1, up=2, strange=3, charm=4, bottom=5, top=6
1639  // The weights are proportional to the square of the
1640  // quark charge (in units of e) and indexed by flavor.
1641  // d=-1/3, u=2/3, s=-1/3, c=2/3, b=-1/3, t=2/3.
1642  // The PDG code for an antiquark is just the negative
1643  // of its partner quark. The charge is also the negative
1644  // of its partner quark charge, so the weights are the
1645  // same for quarks and antiquarks. The PDF are GRSV-standard
1646 
1647  const Double_t weights[] = { -1, 1, 4, 1, 4, 1, 4 };
1648  Double_t polSum = 0;
1649  Double_t unpolSum = 0;
1650 
1651  for (int flavor = 1; flavor <= 6; ++flavor)
1652  {
1653  Double_t polPdf = get_polPDF_NLO(flavor, x, Q2);
1654  polPdf += get_polPDF_NLO(-flavor, x, Q2);
1655  polPdf *= weights[flavor];
1656  polSum += polPdf;
1657  Double_t unpolPdf = get_unpolPDF_NLO(flavor, x, Q2);
1658  unpolPdf += get_unpolPDF_NLO(-flavor, x, Q2);
1659  unpolPdf *= weights[flavor];
1660  unpolSum += unpolPdf;
1661  }
1662 
1663  Double_t A1 = polSum / unpolSum;
1664 
1665  return A1;
1666 
1667 }
1668 
1669 struct PDFWrapper {
1670  typedef double (*PDF)(int,double,double);
1671 
1672  static const PDF pdfs[StPythiaEvent::NPDF];
1673  static PDF pdf;
1674  static int flavor;
1675  static double Q2;
1676 
1677  static void setPdfFlavorQ2(int pdf, int flavor, double Q2)
1678  {
1679  PDFWrapper::pdf = pdfs[pdf];
1680  PDFWrapper::flavor = flavor;
1681  PDFWrapper::Q2 = Q2;
1682  }
1683 
1684  static double eval(double& x) { return pdf(flavor,x,Q2); }
1685 };
1686 
1687 const PDFWrapper::PDF PDFWrapper::pdfs[StPythiaEvent::NPDF] = {
1688  StMCAsymMaker::get_polPDF_LO,
1689  StMCAsymMaker::get_polPDF_NLO, // STD
1690  StMCAsymMaker::get_polPDF_NLO_g0,
1691  StMCAsymMaker::get_polPDF_NLO_gmax,
1692  StMCAsymMaker::get_polPDF_NLO_gmin,
1693  StMCAsymMaker::get_polPDF_NLO_m015,
1694  StMCAsymMaker::get_polPDF_NLO_m030,
1695  StMCAsymMaker::get_polPDF_NLO_m045,
1696  StMCAsymMaker::get_polPDF_NLO_m060,
1697  StMCAsymMaker::get_polPDF_NLO_m075,
1698  StMCAsymMaker::get_polPDF_NLO_m090,
1699  StMCAsymMaker::get_polPDF_NLO_m105,
1700  StMCAsymMaker::get_polPDF_NLO_p030,
1701  StMCAsymMaker::get_polPDF_NLO_p045,
1702  StMCAsymMaker::get_polPDF_NLO_p060,
1703  StMCAsymMaker::get_polPDF_NLO_p070,
1704  StMCAsymMaker::get_polPDF_NLO_GSA,
1705  StMCAsymMaker::get_polPDF_NLO_GSB,
1706  StMCAsymMaker::get_polPDF_NLO_GSC,
1707  StMCAsymMaker::get_polPDF_NLO_DSSV,
1708  StMCAsymMaker::get_polPDF_NLO_LSS1,
1709  StMCAsymMaker::get_polPDF_NLO_LSS2,
1710  StMCAsymMaker::get_polPDF_NLO_LSS3,
1711  StMCAsymMaker::get_polPDF_NLO_AAC1,
1712  StMCAsymMaker::get_polPDF_NLO_AAC2,
1713  StMCAsymMaker::get_polPDF_NLO_AAC3,
1714  StMCAsymMaker::get_polPDF_NLO_BB1,
1715  StMCAsymMaker::get_polPDF_NLO_BB2,
1716  StMCAsymMaker::get_polPDF_NLO_DNS1,
1717  StMCAsymMaker::get_polPDF_NLO_DNS2,
1718  StMCAsymMaker::get_polPDF_NLO_DSSV2009a,
1719  StMCAsymMaker::get_polPDF_NLO_LSS2010_delGpos,
1720  StMCAsymMaker::get_polPDF_NLO_LSS2010_chsign_delG,
1721  StMCAsymMaker::get_polPDF_NLO_BB2010
1722 };
1723 
1724 PDFWrapper::PDF PDFWrapper::pdf = 0;
1725 int PDFWrapper::flavor = 0;
1726 double PDFWrapper::Q2 = 0;
1727 
1728 double StMCAsymMaker::get_polPDF_firstMoment(int pdf, int flavor, double Q2, double xmin, double xmax, double epsilon)
1729 {
1730  PDFWrapper::setPdfFlavorQ2(pdf,flavor,Q2);
1731  return dinteg_(PDFWrapper::eval,xmin,xmax,epsilon);
1732 }
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
StMCAsymMaker::Make
virtual Int_t Make()
Make - this method is called in loop for each event.
Definition: StMCAsymMaker.cxx:202
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
Event
Definition: AgUStep.h:26
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Generated by   doxygen 1.8.5