StRoot: StSpinPool/StMCAsymMaker/StMCAsymMaker.cxx Source File

00001 //*-- Author : Renee Fatemi 
00002 
00003 #include "TFile.h"
00004 #include "StSpinPool/StMCAsymMaker/StMCAsymMaker.h"
00005 #include "StChain.h"
00006 #include "StMcEventMaker/StMcEventMaker.h"
00007 #include "StMcEventTypes.hh"
00008 #include "StMcEvent.hh"
00009 #include "St_DataSet.h"
00010 #include "St_DataSetIter.h"
00011 #include "tables/St_g2t_event_Table.h"
00012 #include "tables/St_particle_Table.h"
00013 #include "tables/St_g2t_pythia_Table.h"
00014 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00015 #include "StMuDSTMaker/COMMON/StMuDst.h"
00016 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00017 
00018 //#include "StMCAsymEvent.h"
00019 #include "StPythiaEvent.h"
00020 
00021 #include "tables/St_particle_Table.h"
00022 
00023 ClassImp(StMCAsymMaker)
00024 
00025 StMCAsymMaker::StMCAsymMaker(const char *name):StMaker(name) {
00026     mEvent = new StPythiaEvent();
00027 }
00028 
00029 StMCAsymMaker::~StMCAsymMaker() {
00030     delete mEvent;
00031 }
00032 
00033 Int_t StMCAsymMaker::Init() {
00034     int iset = 0;
00035     dssvini2009a_(&iset);
00036     return StMaker::Init();
00037 }
00038 
00039 void StMCAsymMaker::Zero() {
00040     pid=-10;
00041     hard_p= -10;
00042     cos_th= -10;
00043     x1= -10;
00044     x2= -10;
00045     s= 0;
00046     t= 0;
00047     u= 0;
00048 
00049     partonic_all=0;
00050     Q2=0;
00051     pid=0;
00052 
00053     df1_NLO_GSA=0;
00054     df2_NLO_GSA=0;
00055     weight_NLO_GSA=0;
00056 
00057     df1_NLO_GSB=0;
00058     df2_NLO_GSB=0;
00059     weight_NLO_GSB=0;
00060 
00061     df1_NLO_GSC=0;
00062     df2_NLO_GSC=0;
00063     weight_NLO_GSC=0;
00064 
00065     df1_LO=0;
00066     df2_LO=0;
00067     f1_LO=0;
00068     f2_LO=0;
00069     weight_LO=0;
00070 
00071     df1_NLO=0;
00072     df2_NLO=0;
00073     f1_NLO=0;
00074     f2_NLO=0;
00075     weight_NLO=0;
00076 
00077     df1_NLO_g0=0;
00078     df2_NLO_g0=0;
00079     weight_NLO_g0=0;
00080 
00081     df1_NLO_gmax=0;
00082     df2_NLO_gmax=0;
00083     weight_NLO_gmax=0;
00084 
00085     df1_NLO_gmin=0;
00086     df2_NLO_gmin=0;
00087     weight_NLO_gmin=0;
00088 
00089     df1_NLO_m015=0;
00090     df2_NLO_m015=0;
00091     weight_NLO_m015=0;
00092 
00093     df1_NLO_m030=0;
00094     df2_NLO_m030=0;
00095     weight_NLO_m030=0;
00096 
00097     df1_NLO_m045=0;
00098     df2_NLO_m045=0;
00099     weight_NLO_m045=0;
00100 
00101     df1_NLO_m060=0;
00102     df2_NLO_m060=0;
00103     weight_NLO_m060=0;
00104 
00105     df1_NLO_m075=0;
00106     df2_NLO_m075=0;
00107     weight_NLO_m075=0;
00108 
00109     df1_NLO_m090=0;
00110     df2_NLO_m090=0;
00111     weight_NLO_m090=0;
00112 
00113     df1_NLO_m105=0;
00114     df2_NLO_m105=0;
00115     weight_NLO_m105=0;
00116 
00117     df1_NLO_p030=0;
00118     df2_NLO_p030=0;
00119     weight_NLO_p030=0;
00120 
00121     df1_NLO_p045=0;
00122     df2_NLO_p045=0;
00123     weight_NLO_p045=0;
00124 
00125     df1_NLO_p060=0;
00126     df2_NLO_p060=0;
00127     weight_NLO_p060=0;
00128 
00129     df1_NLO_p070=0;
00130     df2_NLO_p070=0;
00131     weight_NLO_p070=0;
00132 
00133     df1_NLO_DSSV=0;
00134     df2_NLO_DSSV=0;
00135     weight_NLO_DSSV=0;
00136 
00137     df1_NLO_DSSV2009a=0;
00138     df2_NLO_DSSV2009a=0;
00139     weight_NLO_DSSV2009a=0;
00140 
00141     df1_NLO_LSS1=0;
00142     df2_NLO_LSS1=0;
00143     weight_NLO_LSS1=0;
00144 
00145     df1_NLO_LSS2=0;
00146     df2_NLO_LSS2=0;
00147     weight_NLO_LSS2=0;
00148 
00149     df1_NLO_LSS3=0;
00150     df2_NLO_LSS3=0;
00151     weight_NLO_LSS3=0;
00152 
00153     df1_NLO_AAC1=0;
00154     df2_NLO_AAC1=0;
00155     weight_NLO_AAC1=0;
00156 
00157     df1_NLO_AAC2=0;
00158     df2_NLO_AAC2=0;
00159     weight_NLO_AAC2=0;
00160 
00161     df1_NLO_AAC3=0;
00162     df2_NLO_AAC3=0;
00163     weight_NLO_AAC3=0;
00164 
00165     df1_NLO_BB1=0;
00166     df2_NLO_BB1=0;
00167     weight_NLO_BB1=0;
00168 
00169     df1_NLO_BB2=0;
00170     df2_NLO_BB2=0;
00171     weight_NLO_BB2=0;
00172 
00173     df1_NLO_DNS1=0;
00174     df2_NLO_DNS1=0;
00175     weight_NLO_DNS1=0;
00176 
00177     df1_NLO_DNS2=0;
00178     df2_NLO_DNS2=0;
00179     weight_NLO_DNS2=0;
00180 
00181 }
00182 
00183 void StMCAsymMaker::Clear(const Option_t* c) {
00184     Zero();
00185     mEvent->Clear(c);
00186     StMaker::Clear(c);
00187 }
00188 //_____________________________________________________________________________
00190 Int_t StMCAsymMaker::Make() {
00191     //Get StMcEvent to look at GEANT record
00192     mcEvent = (StMcEvent*)GetDataSet("StMcEvent");
00193     if (!mcEvent) {
00194         LOG_WARN << "No StMcEvent" << endm;
00195         return kStWarn;
00196     }
00197 
00198     //GET EVTID FROM MuDst 
00199     muDstMaker = (StMuDstMaker*)GetMaker("MuDst"); assert(muDstMaker);
00200     StMuDst* dst = muDstMaker->muDst(); assert(dst);
00201     muEvent = dst->event();  assert(muEvent);
00202     StEventInfo &info=muEvent->eventInfo();
00203     evtid=info.id();
00204 
00205     //GET GEANT EVENT
00206     TDataSet *Event = GetDataSet("geant"); //Event->ls(3);
00207 
00208     //GET PYTHIA RECORD from particleTable
00209     TDataSetIter geantDstI(Event);
00210     particleTabPtr = (St_particle  *) geantDstI("particle");
00211 
00212     // if no particle table in the file just skip it
00213     if (particleTabPtr!=0){
00214       
00215       particle_st* particleTable = particleTabPtr->GetTable();//particleTabPtr->Print();
00216       
00217       //GET EVTID and SUBPROCESS ID from struct g2t_event
00218       Pg2t_event=(St_g2t_event *) geantDstI("g2t_event"); //Pg2t_event->Print();
00219       g2t_event_st *g2t_event1=Pg2t_event->GetTable();
00220       geantID= g2t_event1->n_event; 
00221       geantPID= g2t_event1->subprocess_id;
00222       pid=geantPID;
00223       //TEST that geantID==eventID to ensure .geant and .MuDst file are synchronized
00224       assert(evtid==geantID);
00225       
00226       //GET PARTONIC KINEMATICS from struct Pg2t_pythia
00227       Pg2t_pythia=(St_g2t_pythia *) geantDstI("g2t_pythia");// Pg2t_pythia->Print();
00228       g2t_pythia_st *g2t_pythia1=Pg2t_pythia->GetTable();
00229       s= g2t_pythia1-> mand_s;
00230       t= g2t_pythia1-> mand_t;
00231       u= g2t_pythia1-> mand_u;
00232       hard_p= g2t_pythia1->hard_p;
00233       cos_th= g2t_pythia1->cos_th;
00234       x1= g2t_pythia1->bjor_1;
00235       x2= g2t_pythia1->bjor_2;
00236       
00237       //GET FLAVOR AFTER INTIAL RADIATION BEFORE and AFTER SCATTERING
00238       flavor1=particleTable[4].idhep;
00239       flavor2=particleTable[5].idhep;
00240       flavor3=particleTable[6].idhep;
00241       flavor4=particleTable[7].idhep;
00242       
00243       //GET SCATTERED PARTON RECORD
00244       parton1[0]=particleTable[6].idhep;// particle id
00245       parton1[1]=particleTable[6].phep[0];//px
00246       parton1[2]=particleTable[6].phep[1];//py
00247       parton1[3]=particleTable[6].phep[2];//pz
00248       parton1[4]=particleTable[6].phep[3];//E
00249       parton1[5]=particleTable[6].phep[4];//m
00250       parton1[6]=particleTable[6].isthep;//status
00251       parton1[7]=particleTable[6].jmohep[0];//moth1
00252       parton1[8]=particleTable[6].jmohep[1];//moth2
00253       parton1[9]=particleTable[6].jdahep[0];//daughter1
00254       parton1[10]=particleTable[6].jdahep[1];//daughter2
00255       parton2[0]=particleTable[7].idhep;// particle id
00256       parton2[1]=particleTable[7].phep[0];//px
00257       parton2[2]=particleTable[7].phep[1];//py
00258       parton2[3]=particleTable[7].phep[2];//pz
00259       parton2[4]=particleTable[7].phep[3];//E
00260       parton2[5]=particleTable[7].phep[4];//m
00261       parton2[6]=particleTable[7].isthep;//status
00262       parton2[7]=particleTable[7].jmohep[0];//moth1
00263       parton2[8]=particleTable[7].jmohep[1];//moth2
00264       parton2[9]=particleTable[7].jdahep[0];//daughter1
00265       parton2[10]=particleTable[7].jdahep[1];//daughter2
00266       
00267       if (0){//PRINT OUT PYTHIA RECORD
00268         printf("PID/evtid from McEvent = %d,%d; PID/evtid from Table = %d,%d:\n",pid,evtid,geantPID,geantID);
00269         printf("row |   id   |   px   |   py   |   pz   |   E   |   m   | status | moth1 | moth2 | daught1 | daught2 |\n");
00270         for (int i=0; i<particleTabPtr->GetNRows();++i) {
00271           printf("  %d,  %d,  %f,   %f,   %f,   %f,   %f,   %d,   %d,   %d,   %d,   %d\n",i,particleTable[i].idhep, particleTable[i].phep[0],
00272                  particleTable[i].phep[1], particleTable[i].phep[2] , particleTable[i].phep[3], particleTable[i].phep[4], particleTable[i].isthep,
00273                  particleTable[i].jmohep[0], particleTable[i].jmohep[1], particleTable[i].jdahep[0], particleTable[i].jdahep[1]);}
00274         cout << endl;
00275         cout << "flavor1: " << flavor1 << "  flavor2: " << flavor2 << "  flavor3: " << flavor3 << "  flavor4: " << flavor4 << endl << endl;
00276       }
00277       
00278       //Get partonic a_LL, polarized/unpolarized pdfs using Q2 = partonic_pT^2
00279       partonic_all=getPartonicALL(s,t,u,pid,flavor1,flavor2,flavor3,flavor4);
00280       Q2=hard_p*hard_p;
00281       
00282       //NLO GS SCENARIO A
00283       df1_NLO_GSA=get_polPDF_NLO_GSA(flavor1,x1,Q2);
00284       df2_NLO_GSA=get_polPDF_NLO_GSA(flavor2,x2,Q2);
00285       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00286       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00287       weight_NLO_GSA=(df1_NLO_GSA*df2_NLO_GSA*partonic_all)/(f1_NLO*f2_NLO);
00288       
00289       //NLO GS SCENARIO B
00290       df1_NLO_GSB=get_polPDF_NLO_GSB(flavor1,x1,Q2);
00291       df2_NLO_GSB=get_polPDF_NLO_GSB(flavor2,x2,Q2);
00292       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00293       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00294       weight_NLO_GSB=(df1_NLO_GSB*df2_NLO_GSB*partonic_all)/(f1_NLO*f2_NLO);
00295       
00296       //NLO GS SCENARIO C
00297       df1_NLO_GSC=get_polPDF_NLO_GSC(flavor1,x1,Q2);
00298       df2_NLO_GSC=get_polPDF_NLO_GSC(flavor2,x2,Q2);
00299       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00300       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00301       weight_NLO_GSC=(df1_NLO_GSC*df2_NLO_GSC*partonic_all)/(f1_NLO*f2_NLO);
00302       
00303       //LO GRSV
00304       df1_LO=get_polPDF_LO(flavor1,x1,Q2);
00305       df2_LO=get_polPDF_LO(flavor2,x2,Q2);
00306       f1_LO=get_unpolPDF_LO(flavor1,x1,Q2);
00307       f2_LO=get_unpolPDF_LO(flavor2,x2,Q2);
00308       weight_LO=(df1_LO*df2_LO*partonic_all)/(f1_LO*f2_LO);
00309       
00310       //NLO GRSV
00311       df1_NLO=get_polPDF_NLO(flavor1,x1,Q2);
00312       df2_NLO=get_polPDF_NLO(flavor2,x2,Q2);
00313       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00314       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00315       weight_NLO=(df1_NLO*df2_NLO*partonic_all)/(f1_NLO*f2_NLO);
00316       
00317       //NLO_g0 GRSV
00318       df1_NLO_g0=get_polPDF_NLO_g0(flavor1,x1,Q2);
00319       df2_NLO_g0=get_polPDF_NLO_g0(flavor2,x2,Q2);
00320       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00321       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00322       weight_NLO_g0=(df1_NLO_g0*df2_NLO_g0*partonic_all)/(f1_NLO*f2_NLO);
00323       
00324       //NLO_gmax  GRSV
00325       df1_NLO_gmax=get_polPDF_NLO_gmax(flavor1,x1,Q2);
00326       df2_NLO_gmax=get_polPDF_NLO_gmax(flavor2,x2,Q2);
00327       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00328       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00329       weight_NLO_gmax=(df1_NLO_gmax*df2_NLO_gmax*partonic_all)/(f1_NLO*f2_NLO);
00330       
00331       //NLO_gmin  GRSV
00332       df1_NLO_gmin=get_polPDF_NLO_gmin(flavor1,x1,Q2);
00333       df2_NLO_gmin=get_polPDF_NLO_gmin(flavor2,x2,Q2);
00334       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00335       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00336       weight_NLO_gmin=(df1_NLO_gmin*df2_NLO_gmin*partonic_all)/(f1_NLO*f2_NLO);
00337       
00338       //NLO_m015 GRSV
00339       df1_NLO_m015=get_polPDF_NLO_m015(flavor1,x1,Q2);
00340       df2_NLO_m015=get_polPDF_NLO_m015(flavor2,x2,Q2);
00341       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00342       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00343       weight_NLO_m015=(df1_NLO_m015*df2_NLO_m015*partonic_all)/(f1_NLO*f2_NLO);
00344       
00345       //NLO_m030 GRSV
00346       df1_NLO_m030=get_polPDF_NLO_m030(flavor1,x1,Q2);
00347       df2_NLO_m030=get_polPDF_NLO_m030(flavor2,x2,Q2);
00348       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00349       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00350       weight_NLO_m030=(df1_NLO_m030*df2_NLO_m030*partonic_all)/(f1_NLO*f2_NLO);
00351       
00352       //NLO_m045 GRSV
00353       df1_NLO_m045=get_polPDF_NLO_m045(flavor1,x1,Q2);
00354       df2_NLO_m045=get_polPDF_NLO_m045(flavor2,x2,Q2);
00355       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00356       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00357       weight_NLO_m045=(df1_NLO_m045*df2_NLO_m045*partonic_all)/(f1_NLO*f2_NLO);
00358       
00359       //NLO_m060 GRSV
00360       df1_NLO_m060=get_polPDF_NLO_m060(flavor1,x1,Q2);
00361       df2_NLO_m060=get_polPDF_NLO_m060(flavor2,x2,Q2);
00362       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00363       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00364       weight_NLO_m060=(df1_NLO_m060*df2_NLO_m060*partonic_all)/(f1_NLO*f2_NLO);
00365       
00366       //NLO_m075 GRSV
00367       df1_NLO_m075=get_polPDF_NLO_m075(flavor1,x1,Q2);
00368       df2_NLO_m075=get_polPDF_NLO_m075(flavor2,x2,Q2);
00369       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00370       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00371       weight_NLO_m075=(df1_NLO_m075*df2_NLO_m075*partonic_all)/(f1_NLO*f2_NLO);
00372       
00373       //NLO_m090 GRSV
00374       df1_NLO_m090=get_polPDF_NLO_m090(flavor1,x1,Q2);
00375       df2_NLO_m090=get_polPDF_NLO_m090(flavor2,x2,Q2);
00376       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00377       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00378       weight_NLO_m090=(df1_NLO_m090*df2_NLO_m090*partonic_all)/(f1_NLO*f2_NLO);
00379       
00380       //NLO_m105 GRSV
00381       df1_NLO_m105=get_polPDF_NLO_m105(flavor1,x1,Q2);
00382       df2_NLO_m105=get_polPDF_NLO_m105(flavor2,x2,Q2);
00383       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00384       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00385       weight_NLO_m105=(df1_NLO_m105*df2_NLO_m105*partonic_all)/(f1_NLO*f2_NLO);
00386       
00387       //NLO_p030 GRSV
00388       df1_NLO_p030=get_polPDF_NLO_p030(flavor1,x1,Q2);
00389       df2_NLO_p030=get_polPDF_NLO_p030(flavor2,x2,Q2);
00390       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00391       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00392       weight_NLO_p030=(df1_NLO_p030*df2_NLO_p030*partonic_all)/(f1_NLO*f2_NLO);
00393       
00394       //NLO_p045 GRSV
00395       df1_NLO_p045=get_polPDF_NLO_p045(flavor1,x1,Q2);
00396       df2_NLO_p045=get_polPDF_NLO_p045(flavor2,x2,Q2);
00397       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00398       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00399       weight_NLO_p045=(df1_NLO_p045*df2_NLO_p045*partonic_all)/(f1_NLO*f2_NLO);
00400       
00401       //NLO_p060 GRSV
00402       df1_NLO_p060=get_polPDF_NLO_p060(flavor1,x1,Q2);
00403       df2_NLO_p060=get_polPDF_NLO_p060(flavor2,x2,Q2);
00404       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00405       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00406       weight_NLO_p060=(df1_NLO_p060*df2_NLO_p060*partonic_all)/(f1_NLO*f2_NLO);
00407       
00408       //NLO_p070 GRSV
00409       df1_NLO_p070=get_polPDF_NLO_p070(flavor1,x1,Q2);
00410       df2_NLO_p070=get_polPDF_NLO_p070(flavor2,x2,Q2);
00411       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00412       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00413       weight_NLO_p070=(df1_NLO_p070*df2_NLO_p070*partonic_all)/(f1_NLO*f2_NLO);
00414       
00415       //NLO DSSV
00416       df1_NLO_DSSV=get_polPDF_NLO_DSSV(flavor1,x1,Q2);
00417       df2_NLO_DSSV=get_polPDF_NLO_DSSV(flavor2,x2,Q2);
00418       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00419       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00420       weight_NLO_DSSV=(df1_NLO_DSSV*df2_NLO_DSSV*partonic_all)/(f1_NLO*f2_NLO);
00421       
00422       //NLO DSSV 2009a
00423       df1_NLO_DSSV2009a=get_polPDF_NLO_DSSV2009a(flavor1,x1,Q2);
00424       df2_NLO_DSSV2009a=get_polPDF_NLO_DSSV2009a(flavor2,x2,Q2);
00425       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00426       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00427       weight_NLO_DSSV2009a=(df1_NLO_DSSV2009a/f1_NLO)*(df2_NLO_DSSV2009a/f2_NLO)*partonic_all;
00428       
00429       //NLO LSS SCENARIO 1
00430       df1_NLO_LSS1=get_polPDF_NLO_LSS1(flavor1,x1,Q2);
00431       df2_NLO_LSS1=get_polPDF_NLO_LSS1(flavor2,x2,Q2);
00432       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00433       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00434       weight_NLO_LSS1=(df1_NLO_LSS1*df2_NLO_LSS1*partonic_all)/(f1_NLO*f2_NLO);
00435       
00436       //NLO LSS SCENARIO 2
00437       df1_NLO_LSS2=get_polPDF_NLO_LSS2(flavor1,x1,Q2);
00438       df2_NLO_LSS2=get_polPDF_NLO_LSS2(flavor2,x2,Q2);
00439       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00440       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00441       weight_NLO_LSS2=(df1_NLO_LSS2*df2_NLO_LSS2*partonic_all)/(f1_NLO*f2_NLO);
00442       
00443       //NLO LSS SCENARIO 3
00444       df1_NLO_LSS3=get_polPDF_NLO_LSS3(flavor1,x1,Q2);
00445       df2_NLO_LSS3=get_polPDF_NLO_LSS3(flavor2,x2,Q2);
00446       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00447       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00448       weight_NLO_LSS3=(df1_NLO_LSS3*df2_NLO_LSS3*partonic_all)/(f1_NLO*f2_NLO);
00449       
00450       //NLO AAC SCENARIO 1
00451       df1_NLO_AAC1=get_polPDF_NLO_AAC1(flavor1,x1,Q2);
00452       df2_NLO_AAC1=get_polPDF_NLO_AAC1(flavor2,x2,Q2);
00453       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00454       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00455       weight_NLO_AAC1=(df1_NLO_AAC1*df2_NLO_AAC1*partonic_all)/(f1_NLO*f2_NLO);
00456       
00457       //NLO AAC SCENARIO 2
00458       df1_NLO_AAC2=get_polPDF_NLO_AAC2(flavor1,x1,Q2);
00459       df2_NLO_AAC2=get_polPDF_NLO_AAC2(flavor2,x2,Q2);
00460       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00461       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00462       weight_NLO_AAC2=(df1_NLO_AAC2*df2_NLO_AAC2*partonic_all)/(f1_NLO*f2_NLO);
00463       
00464       //NLO AAC SCENARIO 3
00465       df1_NLO_AAC3=get_polPDF_NLO_AAC3(flavor1,x1,Q2);
00466       df2_NLO_AAC3=get_polPDF_NLO_AAC3(flavor2,x2,Q2);
00467       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00468       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00469       weight_NLO_AAC3=(df1_NLO_AAC3*df2_NLO_AAC3*partonic_all)/(f1_NLO*f2_NLO);
00470       
00471       //NLO BB SCENARIO 1
00472       df1_NLO_BB1=get_polPDF_NLO_BB1(flavor1,x1,Q2);
00473       df2_NLO_BB1=get_polPDF_NLO_BB1(flavor2,x2,Q2);
00474       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00475       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00476       weight_NLO_BB1=(df1_NLO_BB1*df2_NLO_BB1*partonic_all)/(f1_NLO*f2_NLO);
00477       
00478       //NLO BB SCENARIO 2
00479       df1_NLO_BB2=get_polPDF_NLO_BB2(flavor1,x1,Q2);
00480       df2_NLO_BB2=get_polPDF_NLO_BB2(flavor2,x2,Q2);
00481       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00482       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00483       weight_NLO_BB2=(df1_NLO_BB2*df2_NLO_BB2*partonic_all)/(f1_NLO*f2_NLO);
00484       
00485       //NLO DNS SCENARIO 1
00486       df1_NLO_DNS1=get_polPDF_NLO_DNS1(flavor1,x1,Q2);
00487       df2_NLO_DNS1=get_polPDF_NLO_DNS1(flavor2,x2,Q2);
00488       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00489       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00490       weight_NLO_DNS1=(df1_NLO_DNS1*df2_NLO_DNS1*partonic_all)/(f1_NLO*f2_NLO);
00491       
00492       //NLO DNS SCENARIO 2
00493       df1_NLO_DNS2=get_polPDF_NLO_DNS2(flavor1,x1,Q2);
00494       df2_NLO_DNS2=get_polPDF_NLO_DNS2(flavor2,x2,Q2);
00495       f1_NLO=get_unpolPDF_NLO(flavor1,x1,Q2);
00496       f2_NLO=get_unpolPDF_NLO(flavor2,x2,Q2);
00497       weight_NLO_DNS2=(df1_NLO_DNS2*df2_NLO_DNS2*partonic_all)/(f1_NLO*f2_NLO);
00498       
00499       if (0) {
00500         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);
00501         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);
00502         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);
00503         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);
00504         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);
00505         
00506         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);
00507         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);
00508         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);
00509         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);
00510         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);
00511         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);
00512         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);
00513         
00514         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);
00515         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);
00516         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);
00517         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);
00518         
00519         
00520         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);
00521         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);
00522         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);
00523         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);
00524         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);
00525         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);
00526         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);
00527         
00528         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);
00529         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);
00530         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);
00531         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);
00532       }
00533     
00534  
00535     fillPythiaEvent(mEvent);
00536 
00537     }
00538 
00539     return kStOK;
00540 }
00541 
00542 void StMCAsymMaker::fillPythiaEvent(StPythiaEvent* pythia)
00543 {
00544     g2t_event_st* eventTable = Pg2t_event->GetTable();
00545 
00546     pythia->setRunId(eventTable->n_run);
00547     pythia->setEventId(eventTable->n_event);
00548 
00549     if (mcEvent && mcEvent->primaryVertex())
00550         pythia->setVertex(mcEvent->primaryVertex()->position().xyz());
00551 
00552     g2t_pythia_st* pythiaTable = Pg2t_pythia->GetTable();
00553 
00554     pythia->setProcessId(pythiaTable->subprocess_id);
00555     pythia->setS(pythiaTable->mand_s);
00556     pythia->setT(pythiaTable->mand_t);
00557     pythia->setU(pythiaTable->mand_u);
00558     pythia->setPt(pythiaTable->hard_p);
00559     pythia->setCosTheta(pythiaTable->cos_th);
00560     pythia->setX1(pythiaTable->bjor_1);
00561     pythia->setX2(pythiaTable->bjor_2);
00562     pythia->setMstu72(pythiaTable->mstu72);
00563     pythia->setMstu73(pythiaTable->mstu73);
00564     pythia->setMstp111(pythiaTable->mstp111);
00565     pythia->setPartonALL(partonic_all);
00566     
00567     pythia->setDF1(StPythiaEvent::LO, df1_LO);
00568     pythia->setDF1(StPythiaEvent::NLO, df1_NLO);
00569     pythia->setDF1(StPythiaEvent::ZERO, df1_NLO_g0);
00570     pythia->setDF1(StPythiaEvent::MAX, df1_NLO_gmax);
00571     pythia->setDF1(StPythiaEvent::MIN, df1_NLO_gmin);
00572     pythia->setDF1(StPythiaEvent::M015, df1_NLO_m015);
00573     pythia->setDF1(StPythiaEvent::M030, df1_NLO_m030);
00574     pythia->setDF1(StPythiaEvent::M045, df1_NLO_m045);
00575     pythia->setDF1(StPythiaEvent::M060, df1_NLO_m060);
00576     pythia->setDF1(StPythiaEvent::M075, df1_NLO_m075);
00577     pythia->setDF1(StPythiaEvent::M090, df1_NLO_m090);
00578     pythia->setDF1(StPythiaEvent::M105, df1_NLO_m105);
00579     pythia->setDF1(StPythiaEvent::P030, df1_NLO_p030);
00580     pythia->setDF1(StPythiaEvent::P045, df1_NLO_p045);
00581     pythia->setDF1(StPythiaEvent::P060, df1_NLO_p060);
00582     pythia->setDF1(StPythiaEvent::P070, df1_NLO_p070);
00583     pythia->setDF1(StPythiaEvent::GS_NLOA, df1_NLO_GSA);
00584     pythia->setDF1(StPythiaEvent::GS_NLOB, df1_NLO_GSB);
00585     pythia->setDF1(StPythiaEvent::GS_NLOC, df1_NLO_GSC);
00586     pythia->setDF1(StPythiaEvent::DSSV, df1_NLO_DSSV);
00587     pythia->setDF1(StPythiaEvent::DSSV2009a, df1_NLO_DSSV2009a);
00588     pythia->setDF1(StPythiaEvent::LSS1, df1_NLO_LSS1);
00589     pythia->setDF1(StPythiaEvent::LSS2, df1_NLO_LSS2);
00590     pythia->setDF1(StPythiaEvent::LSS3, df1_NLO_LSS3);
00591     pythia->setDF1(StPythiaEvent::AAC1, df1_NLO_AAC1);
00592     pythia->setDF1(StPythiaEvent::AAC2, df1_NLO_AAC2);
00593     pythia->setDF1(StPythiaEvent::AAC3, df1_NLO_AAC3);
00594     pythia->setDF1(StPythiaEvent::BB1, df1_NLO_BB1);
00595     pythia->setDF1(StPythiaEvent::BB2, df1_NLO_BB2);
00596     pythia->setDF1(StPythiaEvent::DNS1, df1_NLO_DNS1);
00597     pythia->setDF1(StPythiaEvent::DNS2, df1_NLO_DNS2);
00598 
00599     pythia->setDF2(StPythiaEvent::LO, df2_LO);
00600     pythia->setDF2(StPythiaEvent::NLO, df2_NLO);
00601     pythia->setDF2(StPythiaEvent::ZERO, df2_NLO_g0);
00602     pythia->setDF2(StPythiaEvent::MAX, df2_NLO_gmax);
00603     pythia->setDF2(StPythiaEvent::MIN, df2_NLO_gmin);
00604     pythia->setDF2(StPythiaEvent::M015, df2_NLO_m015);
00605     pythia->setDF2(StPythiaEvent::M030, df2_NLO_m030);
00606     pythia->setDF2(StPythiaEvent::M045, df2_NLO_m045);
00607     pythia->setDF2(StPythiaEvent::M060, df2_NLO_m060);
00608     pythia->setDF2(StPythiaEvent::M075, df2_NLO_m075);
00609     pythia->setDF2(StPythiaEvent::M090, df2_NLO_m090);
00610     pythia->setDF2(StPythiaEvent::M105, df2_NLO_m105);
00611     pythia->setDF2(StPythiaEvent::P030, df2_NLO_p030);
00612     pythia->setDF2(StPythiaEvent::P045, df2_NLO_p045);
00613     pythia->setDF2(StPythiaEvent::P060, df2_NLO_p060);
00614     pythia->setDF2(StPythiaEvent::P070, df2_NLO_p070);
00615     pythia->setDF2(StPythiaEvent::GS_NLOA, df2_NLO_GSA);
00616     pythia->setDF2(StPythiaEvent::GS_NLOB, df2_NLO_GSB);
00617     pythia->setDF2(StPythiaEvent::GS_NLOC, df2_NLO_GSC);
00618     pythia->setDF2(StPythiaEvent::DSSV, df2_NLO_DSSV);
00619     pythia->setDF2(StPythiaEvent::DSSV2009a, df2_NLO_DSSV2009a);
00620     pythia->setDF2(StPythiaEvent::LSS1, df2_NLO_LSS1);
00621     pythia->setDF2(StPythiaEvent::LSS2, df2_NLO_LSS2);
00622     pythia->setDF2(StPythiaEvent::LSS3, df2_NLO_LSS3);
00623     pythia->setDF2(StPythiaEvent::AAC1, df2_NLO_AAC1);
00624     pythia->setDF2(StPythiaEvent::AAC2, df2_NLO_AAC2);
00625     pythia->setDF2(StPythiaEvent::AAC3, df2_NLO_AAC3);
00626     pythia->setDF2(StPythiaEvent::BB1, df2_NLO_BB1);
00627     pythia->setDF2(StPythiaEvent::BB2, df2_NLO_BB2);
00628     pythia->setDF2(StPythiaEvent::DNS1, df2_NLO_DNS1);
00629     pythia->setDF2(StPythiaEvent::DNS2, df2_NLO_DNS2);
00630     
00631     pythia->setF1(StPythiaEvent::LO, f1_LO);
00632     pythia->setF1(StPythiaEvent::NLO, f1_NLO);
00633     
00634     pythia->setF2(StPythiaEvent::LO, f2_LO);
00635     pythia->setF2(StPythiaEvent::NLO, f2_NLO);
00636 
00637     particle_st* particleTable = particleTabPtr->GetTable();
00638 
00639     for (int i = 4; i < 8; ++i)
00640       pythia->addParticle(TParticle(particleTable[i].idhep,
00641                                     particleTable[i].isthep,
00642                                     particleTable[i].jmohep[0],
00643                                     particleTable[i].jmohep[1],
00644                                     particleTable[i].jdahep[0],
00645                                     particleTable[i].jdahep[1],
00646                                     TLorentzVector(particleTable[i].phep),
00647                                     TLorentzVector(particleTable[i].vhep)));
00648 
00649 }
00650 
00651 //Gehrmann-Stirling NL0 Set A
00652 Double_t StMCAsymMaker::get_polPDF_NLO_GSA(int flavor, double x, double Q2){
00653 
00654     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00655     double pdf=1000;
00656     int polset_NLO_GSA=201;
00657     int polid=0;
00658 
00659     //    cout<<"get_polPDF_NLO_GSA: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;   
00660     if ((Q2>=1.0)&&(Q2<=1E+06))  polar_(&polset_NLO_GSA, &x, &Q2, parpol, &polid);
00661     //    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;
00662 
00663   
00664     //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
00665     if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
00666     //parpol[1] = dv - dbar
00667     if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
00668     //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
00669     if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark 
00670     //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
00671     if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
00672     //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
00673     if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
00674     if (flavor==21) pdf=parpol[5];     //gluon
00675     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
00676 
00677   return pdf;
00678 }
00679 
00680 //Gehrmann-Stirling NL0 Set B
00681 Double_t StMCAsymMaker::get_polPDF_NLO_GSB(int flavor, double x, double Q2){
00682 
00683     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00684     double pdf=1000;
00685     int polset_NLO_GSB=202;
00686     int polid=0;
00687 
00688     //cout<<"get_polPDF_NLO_GSB: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00689     if ((Q2>=1.0)&&(Q2<=1E+06))  polar_(&polset_NLO_GSB, &x, &Q2, parpol, &polid);
00690     //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;
00691     
00692   
00693     //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
00694     if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
00695     //parpol[1] = dv - dbar
00696     if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
00697     //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
00698     if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark 
00699     //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
00700     if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
00701     //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
00702     if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
00703     if (flavor==21) pdf=parpol[5];     //gluon
00704     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
00705 
00706   return pdf;
00707 }
00708 
00709 //Gehrmann-Stirling NL0 Set C
00710 Double_t StMCAsymMaker::get_polPDF_NLO_GSC(int flavor, double x, double Q2){
00711 
00712     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00713     double pdf=1000;
00714     int polset_NLO_GSC=203;
00715     int polid=0;
00716 
00717     //cout<<"get_polPDF_NLO_GSC: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00718     if ((Q2>=1.0)&&(Q2<=1E+06))  polar_(&polset_NLO_GSC, &x, &Q2, parpol, &polid);
00719     //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;
00720    
00721     //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
00722     if (flavor==2) pdf=parpol[0]+parpol[2]; //uv + usea
00723     //parpol[1] = dv - dbar
00724     if (flavor==1) pdf=parpol[1]+parpol[3]; //dv + dsea quark
00725     //parpol[2] = 1/2 usea= 1/2(usea +ubar) = ubar
00726     if (flavor==-2) pdf=parpol[2];//ubar==1/2usea quark 
00727     //parpol[3] = 1/2 dsea= 1/2(dsea +dbar) = dbar
00728     if (flavor==-1) pdf=parpol[3];//dbar==1/2dsea quark
00729     //parpol[4] = 1/2 strangesea = 1/2*(ssea+sbar) = sbar
00730     if (abs(flavor)==3) pdf=parpol[4]; //s==1/2sbar quark
00731     if (flavor==21) pdf=parpol[5];     //gluon
00732     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];//c+b+t == 1/2 strangesea
00733 
00734 
00735   return pdf;
00736 }
00737 
00738 
00739 
00740 
00741 
00742 //GRSV LO standard
00743 Double_t StMCAsymMaker::get_polPDF_LO(int flavor, double x, double Q2){
00744 
00745     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00746     double pdf=1000;
00747     int polset_LO=101;
00748     int polid=0;
00749     //cout<<"get_polPDF_LO: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00750 
00751     if (Q2>=0.8&&Q2<=1.0e6)  polar_(&polset_LO, &x, &Q2, parpol, &polid);
00752     //cout <<"getpolPDF_LO:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00753 
00754     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00755     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00756     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00757     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00758     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00759     if (flavor==21) pdf=parpol[5];     //gluon
00760     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00761 
00762     return pdf;
00763 }
00764 
00765 
00766 //GRSV NLO standard
00767 Double_t StMCAsymMaker::get_polPDF_NLO(int flavor, double x, double Q2){
00768 
00769     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00770     double pdf=1000;
00771     int polset_NLO=102;
00772     int polid=0;
00773     //cout<<"get_polPDF_NLO: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00774 
00775     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO, &x, &Q2, parpol, &polid);
00776     //cout <<"get_polPDF_NLO:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00777 
00778     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00779     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00780     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00781     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00782     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00783     if (flavor==21) pdf=parpol[5];     //gluon
00784     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00785 
00786     return pdf;
00787 }
00788 
00789 //GRSV NLO g0
00790 Double_t StMCAsymMaker::get_polPDF_NLO_g0(int flavor, double x, double Q2){
00791 
00792     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00793     double pdf=1000;
00794     int polset_NLO_g0=103;
00795     int polid=0;
00796     //cout<<"get_polPDF_NLO_g0: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00797 
00798     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_g0, &x, &Q2, parpol, &polid);
00799     //cout <<"getpolPDF_NLO_g0:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00800 
00801     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00802     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00803     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00804     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00805     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00806     if (flavor==21) pdf=parpol[5];     //gluon
00807     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00808 
00809     return pdf;
00810 }
00811 
00812 
00813 //GRSV NLO gmax
00814 Double_t StMCAsymMaker::get_polPDF_NLO_gmax(int flavor, double x, double Q2){
00815 
00816     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00817     double pdf=1000;
00818     int polset_NLO_gmax=104;
00819     int polid=0;
00820     //cout<<"get_polPDF_NLO_gmax: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00821 
00822     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_gmax, &x, &Q2, parpol, &polid);
00823     //cout <<"getpolPDF_NLO_gmax:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00824 
00825     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00826     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00827     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00828     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00829     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00830     if (flavor==21) pdf=parpol[5];     //gluon
00831     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00832 
00833     return pdf;
00834 }
00835 
00836 
00837 //GRSV NLO gmin
00838 Double_t StMCAsymMaker::get_polPDF_NLO_gmin(int flavor, double x, double Q2){
00839 
00840     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00841     double pdf=1000;
00842     int polset_NLO_gmin=105;
00843     int polid=0;
00844     //cout<<"get_polPDF_NLO_gmin: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00845 
00846     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_gmin, &x, &Q2, parpol, &polid);
00847     //cout <<"getpolPDF_NLO_gmin:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00848 
00849     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00850     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00851     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00852     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00853     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00854     if (flavor==21) pdf=parpol[5];     //gluon
00855     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00856 
00857     return pdf;
00858 }
00859 
00860 //GRSV NLO m015
00861 Double_t StMCAsymMaker::get_polPDF_NLO_m015(int flavor, double x, double Q2){
00862 
00863     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00864     double pdf=1000;
00865     int polset_NLO_m015=106;
00866     int polid=0;
00867     //cout<<"get_polPDF_NLO_m015: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00868 
00869     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m015, &x, &Q2, parpol, &polid);
00870     //cout <<"getpolPDF_NLO_m015:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00871 
00872     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00873     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00874     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00875     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00876     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00877     if (flavor==21) pdf=parpol[5];     //gluon
00878     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00879 
00880     return pdf;
00881 }
00882 
00883 //GRSV NLO m030
00884 Double_t StMCAsymMaker::get_polPDF_NLO_m030(int flavor, double x, double Q2){
00885 
00886     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00887     double pdf=1000;
00888     int polset_NLO_m030=107;
00889     int polid=0;
00890     //cout<<"get_polPDF_NLO_m030: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00891 
00892     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m030, &x, &Q2, parpol, &polid);
00893     //cout <<"getpolPDF_NLO_m030:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00894 
00895     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00896     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00897     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00898     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00899     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00900     if (flavor==21) pdf=parpol[5];     //gluon
00901     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00902 
00903     return pdf;
00904 }
00905 
00906 //GRSV NLO m045
00907 Double_t StMCAsymMaker::get_polPDF_NLO_m045(int flavor, double x, double Q2){
00908 
00909     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00910     double pdf=1000;
00911     int polset_NLO_m045=108;
00912     int polid=0;
00913     //cout<<"get_polPDF_NLO_m045: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00914 
00915     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m045, &x, &Q2, parpol, &polid);
00916     //cout <<"getpolPDF_NLO_m045:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00917 
00918     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00919     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00920     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00921     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00922     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00923     if (flavor==21) pdf=parpol[5];     //gluon
00924     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00925 
00926     return pdf;
00927 }
00928 
00929 //GRSV NLO m060
00930 Double_t StMCAsymMaker::get_polPDF_NLO_m060(int flavor, double x, double Q2){
00931 
00932     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00933     double pdf=1000;
00934     int polset_NLO_m060=109;
00935     int polid=0;
00936     //cout<<"get_polPDF_NLO_m060: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00937 
00938     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m060, &x, &Q2, parpol, &polid);
00939     //cout <<"getpolPDF_NLO_m060:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00940 
00941     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00942     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00943     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00944     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00945     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00946     if (flavor==21) pdf=parpol[5];     //gluon
00947     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00948 
00949     return pdf;
00950 }
00951 
00952 //GRSV NLO m075
00953 Double_t StMCAsymMaker::get_polPDF_NLO_m075(int flavor, double x, double Q2){
00954 
00955     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00956     double pdf=1000;
00957     int polset_NLO_m075=110;
00958     int polid=0;
00959     //cout<<"get_polPDF_NLO_m075: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00960 
00961     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m075, &x, &Q2, parpol, &polid);
00962     //cout <<"getpolPDF_NLO_m075:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00963 
00964     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00965     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00966     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00967     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00968     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00969     if (flavor==21) pdf=parpol[5];     //gluon
00970     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00971 
00972     return pdf;
00973 }
00974 
00975 //GRSV NLO m090
00976 Double_t StMCAsymMaker::get_polPDF_NLO_m090(int flavor, double x, double Q2){
00977 
00978     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
00979     double pdf=1000;
00980     int polset_NLO_m090=111;
00981     int polid=0;
00982     //cout<<"get_polPDF_NLO_m090: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
00983 
00984     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m090, &x, &Q2, parpol, &polid);
00985     //cout <<"getpolPDF_NLO_m090:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
00986 
00987     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
00988     if (flavor==2) pdf=parpol[0];      //uv + usea quark
00989     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
00990     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
00991     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
00992     if (flavor==21) pdf=parpol[5];     //gluon
00993     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
00994 
00995     return pdf;
00996 }
00997 
00998 //GRSV NLO m105
00999 Double_t StMCAsymMaker::get_polPDF_NLO_m105(int flavor, double x, double Q2){
01000 
01001     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01002     double pdf=1000;
01003     int polset_NLO_m105=112;
01004     int polid=0;
01005     //cout<<"get_polPDF_NLO_m105: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01006 
01007     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_m105, &x, &Q2, parpol, &polid);
01008     //cout <<"getpolPDF_NLO_m105:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01009 
01010     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01011     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01012     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01013     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01014     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01015     if (flavor==21) pdf=parpol[5];     //gluon
01016     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01017 
01018     return pdf;
01019 }
01020 
01021 //GRSV NLO p030
01022 Double_t StMCAsymMaker::get_polPDF_NLO_p030(int flavor, double x, double Q2){
01023 
01024     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01025     double pdf=1000;
01026     int polset_NLO_p030=113;
01027     int polid=0;
01028     //cout<<"get_polPDF_NLO_p030: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01029 
01030     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p030, &x, &Q2, parpol, &polid);
01031     //cout <<"getpolPDF_NLO_p030:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01032 
01033     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01034     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01035     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01036     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01037     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01038     if (flavor==21) pdf=parpol[5];     //gluon
01039     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01040 
01041     return pdf;
01042 }
01043 
01044 //GRSV NLO p045
01045 Double_t StMCAsymMaker::get_polPDF_NLO_p045(int flavor, double x, double Q2){
01046 
01047     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01048     double pdf=1000;
01049     int polset_NLO_p045=114;
01050     int polid=0;
01051     //cout<<"get_polPDF_NLO_p045: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01052 
01053     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p045, &x, &Q2, parpol, &polid);
01054     //cout <<"getpolPDF_NLO_p045:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01055 
01056     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01057     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01058     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01059     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01060     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01061     if (flavor==21) pdf=parpol[5];     //gluon
01062     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01063 
01064     return pdf;
01065 }
01066 
01067 //GRSV NLO p060
01068 Double_t StMCAsymMaker::get_polPDF_NLO_p060(int flavor, double x, double Q2){
01069 
01070     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01071     double pdf=1000;
01072     int polset_NLO_p060=115;
01073     int polid=0;
01074     //cout<<"get_polPDF_NLO_p060: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01075 
01076     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p060, &x, &Q2, parpol, &polid);
01077     //cout <<"getpolPDF_NLO_p060:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01078 
01079     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01080     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01081     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01082     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01083     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01084     if (flavor==21) pdf=parpol[5];     //gluon
01085     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01086 
01087     return pdf;
01088 }
01089 
01090 //GRSV NLO p070
01091 Double_t StMCAsymMaker::get_polPDF_NLO_p070(int flavor, double x, double Q2){
01092 
01093     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01094     double pdf=1000;
01095     int polset_NLO_p070=116;
01096     int polid=0;
01097     //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01098 
01099     if (Q2>=0.8&&Q2<=1.0e6) polar_(&polset_NLO_p070, &x, &Q2, parpol, &polid);
01100     //cout <<"getpolPDF_NLO_p070:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01101 
01102     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01103     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01104     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01105     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01106     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01107     if (flavor==21) pdf=parpol[5];     //gluon
01108     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01109 
01110     return pdf;
01111 }
01112 
01113 //DSSV NLO
01114 Double_t StMCAsymMaker::get_polPDF_NLO_DSSV(int flavor, double x, double Q2){
01115 
01116     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01117     double pdf=1000;
01118     int polset_NLO_DSSV=301;
01119     int polid=0;
01120     //cout<<"get_polPDF_NLO_DSSV: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01121 
01122     if (Q2>=1.0&&Q2<=1.0e5) polar_(&polset_NLO_DSSV, &x, &Q2, parpol, &polid);
01123     //cout <<"getpolPDF_NLO_DSSV:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01124 
01125     if (flavor==1) pdf=parpol[1]+parpol[2];      //dv + dsea quark
01126     if (flavor==2) pdf=parpol[0]+parpol[3];      //uv + usea quark
01127     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01128     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01129     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01130     if (flavor==21) pdf=parpol[5];     //gluon
01131     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01132 
01133     return pdf;
01134 }
01135 
01136 //DSSV 2009a NLO
01137 Double_t StMCAsymMaker::get_polPDF_NLO_DSSV2009a(int flavor, double x, double Q2){
01138 
01139     double duv=0,ddv=0,dubar=0,ddbar=0,dstr=0,dglu=0;
01140 
01141     if (x>=1.0e-5&&x<=1.0&&Q2>=1.0&&Q2<=1.0e5) {
01142       dssvfit2009a_(&x,&Q2,&duv,&ddv,&dubar,&ddbar,&dstr,&dglu);
01143       switch (flavor) {
01144       case  1: return (ddv+ddbar)/x; //dv + dsea quark
01145       case  2: return (duv+dubar)/x; //uv + usea quark
01146       case -1: return ddbar/x;  //dbar==dsea quark
01147       case -2: return dubar/x;  //ubar==usea quark
01148       case  3:
01149       case -3: return dstr/x;   //s==sbar quark
01150       case 21: return dglu/x;   //gluon
01151       case  4:
01152       case -4:
01153       case  5:
01154       case -5:
01155       case  6:
01156       case -6: return dstr/x;
01157       }
01158     }
01159     return 1000;
01160 }
01161 
01162 //LSS NLO Scenario 1
01163 Double_t StMCAsymMaker::get_polPDF_NLO_LSS1(int flavor, double x, double Q2){
01164 
01165     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01166     double pdf=1000;
01167     int polset_NLO_LSS1=401;
01168     int polid=0;
01169     //    cout<<"get_polPDF_NLO_LSS1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01170 
01171     if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS1, &x, &Q2, parpol, &polid);
01172     //    cout <<"getpolPDF_NLO_LSS1:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01173 
01174     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01175     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01176     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01177     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01178     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01179     if (flavor==21) pdf=parpol[5];     //gluon
01180     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01181 
01182     return pdf;
01183 }
01184 
01185 //LSS NLO Scenario 2
01186 Double_t StMCAsymMaker::get_polPDF_NLO_LSS2(int flavor, double x, double Q2){
01187 
01188     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01189     double pdf=1000;
01190     int polset_NLO_LSS2=402;
01191     int polid=0;
01192     //    cout<<"get_polPDF_NLO_LSS2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01193 
01194     if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS2, &x, &Q2, parpol, &polid);
01195     //    cout <<"getpolPDF_NLO_LSS2:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01196 
01197     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01198     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01199     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01200     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01201     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01202     if (flavor==21) pdf=parpol[5];     //gluon
01203     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01204 
01205     return pdf;
01206 }
01207 
01208 //LSS NLO Scenario 3
01209 Double_t StMCAsymMaker::get_polPDF_NLO_LSS3(int flavor, double x, double Q2){
01210 
01211     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01212     double pdf=1000;
01213     int polset_NLO_LSS3=403;
01214     int polid=0;
01215     //    cout<<"get_polPDF_NLO_LSS3: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01216 
01217     if (Q2>=1.0&&Q2<=5.8e5) polar_(&polset_NLO_LSS3, &x, &Q2, parpol, &polid);
01218     //    cout <<"getpolPDF_NLO_LSS3:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01219 
01220     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01221     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01222     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01223     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01224     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01225     if (flavor==21) pdf=parpol[5];     //gluon
01226     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01227 
01228     return pdf;
01229 }
01230 
01231 //AAC NLO Scenario 1
01232 Double_t StMCAsymMaker::get_polPDF_NLO_AAC1(int flavor, double x, double Q2){
01233 
01234     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01235     double pdf=1000;
01236     int polset_NLO_AAC1=501;
01237     int polid=0;
01238     //cout<<"get_polPDF_NLO_AAC1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01239 
01240     if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC1, &x, &Q2, parpol, &polid);
01241     //cout <<"getpolPDF_NLO_AAC1:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01242 
01243     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01244     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01245     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01246     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01247     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01248     if (flavor==21) pdf=parpol[5];     //gluon
01249     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01250 
01251     return pdf;
01252 }
01253 
01254 //AAC NLO Scenario 2
01255 Double_t StMCAsymMaker::get_polPDF_NLO_AAC2(int flavor, double x, double Q2){
01256 
01257     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01258     double pdf=1000;
01259     int polset_NLO_AAC2=502;
01260     int polid=0;
01261     //cout<<"get_polPDF_NLO_AAC2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01262 
01263     if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC2, &x, &Q2, parpol, &polid);
01264     //cout <<"getpolPDF_NLO_AAC2:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01265 
01266     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01267     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01268     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01269     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01270     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01271     if (flavor==21) pdf=parpol[5];     //gluon
01272     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01273 
01274     return pdf;
01275 }
01276 
01277 //AAC NLO Scenario 3
01278 Double_t StMCAsymMaker::get_polPDF_NLO_AAC3(int flavor, double x, double Q2){
01279 
01280     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01281     double pdf=1000;
01282     int polset_NLO_AAC3=503;
01283     int polid=0;
01284     //cout<<"get_polPDF_NLO_AAC3: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01285 
01286     if (Q2>=1.0&&Q2<=1.0e8) polar_(&polset_NLO_AAC3, &x, &Q2, parpol, &polid);
01287     //cout <<"getpolPDF_NLO_AAC3:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01288 
01289     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01290     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01291     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01292     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01293     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01294     if (flavor==21) pdf=parpol[5];     //gluon
01295     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01296 
01297     return pdf;
01298 }
01299 
01300 //BB NLO Scenario 1
01301 Double_t StMCAsymMaker::get_polPDF_NLO_BB1(int flavor, double x, double Q2){
01302 
01303     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01304     double pdf=1000;
01305     int polset_NLO_BB1=601;
01306     int polid=0;
01307     //cout<<"get_polPDF_NLO_BB1: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01308 
01309     if (Q2>=1.0&&Q2<=1.0e6) polar_(&polset_NLO_BB1, &x, &Q2, parpol, &polid);
01310     //cout <<"get_polPDF_NLO_BB1:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01311 
01312     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01313     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01314     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01315     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01316     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01317     if (flavor==21) pdf=parpol[5];     //gluon
01318     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01319 
01320     return pdf;
01321 }
01322 
01323 //BB NLO Scenario 2
01324 Double_t StMCAsymMaker::get_polPDF_NLO_BB2(int flavor, double x, double Q2){
01325 
01326     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01327     double pdf=1000;
01328     int polset_NLO_BB2=602;
01329     int polid=0;
01330     //cout<<"get_polPDF_NLO_BB2: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01331 
01332     if (Q2>=1.0&&Q2<=1.0e6) polar_(&polset_NLO_BB2, &x, &Q2, parpol, &polid);
01333     //cout <<"get_polPDF_NLO_BB2:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01334 
01335     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01336     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01337     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01338     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01339     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01340     if (flavor==21) pdf=parpol[5];     //gluon
01341     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01342 
01343     return pdf;
01344 }
01345 
01346 //DNS NLO Scenario 1
01347 Double_t StMCAsymMaker::get_polPDF_NLO_DNS1(int flavor, double x, double Q2){
01348 
01349     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01350     double pdf=1000;
01351     int polset_NLO_DNS1=701;
01352     int polid=0;
01353     //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01354 
01355     if (Q2>=1.0&&Q2<=5.0e4) polar_(&polset_NLO_DNS1, &x, &Q2, parpol, &polid);
01356     //cout <<"getpolPDF_NLO_p070:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01357 
01358     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01359     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01360     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01361     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01362     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01363     if (flavor==21) pdf=parpol[5];     //gluon
01364     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01365 
01366     return pdf;
01367 }
01368 
01369 //DNS NLO Scenario 2
01370 Double_t StMCAsymMaker::get_polPDF_NLO_DNS2(int flavor, double x, double Q2){
01371 
01372     double parpol[6]={0.0,0.0,0.0,0.0,0.0,0.0};
01373     double pdf=1000;
01374     int polset_NLO_DNS2=702;
01375     int polid=0;
01376     //cout<<"get_polPDF_NLO_p070: flavor="<<flavor<<" x="<<x<<" Q2="<<Q2<<" id="<<polid<<endl;
01377 
01378     if (Q2>=1.0&&Q2<=5.0e4) polar_(&polset_NLO_DNS2, &x, &Q2, parpol, &polid);
01379     //cout <<"getpolPDF_NLO_p070:  U="<<parpol[0]<<" D="<<parpol[1]<<" UB="<<parpol[2]<<" DB="<<parpol[3]<<" ST="<<parpol[4]<<" GL="<<parpol[5]<<endl;
01380 
01381     if (flavor==1) pdf=parpol[1];      //dv + dsea quark
01382     if (flavor==2) pdf=parpol[0];      //uv + usea quark
01383     if (flavor==-1) pdf=parpol[3];     //dbar==dsea quark
01384     if (flavor==-2) pdf=parpol[2];     //ubar==usea quark
01385     if (abs(flavor)==3) pdf=parpol[4]; //s==sbar quark
01386     if (flavor==21) pdf=parpol[5];     //gluon
01387     if ((abs(flavor)>=4)&&(abs(flavor)<=6)) pdf=parpol[4];
01388 
01389     return pdf;
01390 }
01391 
01392 //Returns unpolarized CTEQ
01393 Double_t StMCAsymMaker::get_unpolPDF_LO(int flavor, double x, double Q2){
01394 
01395     Double_t pdf=0.0;
01396     Int_t iset=3;//LO
01397     Int_t er=0;
01398     Int_t fl=10;
01399 
01400     if (flavor==1) fl=2;
01401     if (flavor==2) fl=1;
01402     if (flavor==-1) fl=-2;
01403     if (flavor==-2) fl=-1;
01404     if (flavor==21) fl=0;
01405     if (flavor==3) fl=3;
01406     if (flavor==-3) fl=-3;
01407     if (flavor==4) fl=4;
01408     if (flavor==-4) fl=-4;
01409     if (flavor==5) fl=5;
01410     if (flavor==-5) fl=-5;
01411 
01412     double Q=pow(Q2,0.5);
01413     pdf=ctq5pd_(&iset,&fl,&x,&Q,&er);
01414     if (er!=0) pdf=0.0;
01415 
01416     return pdf;
01417 }
01418 
01419 Double_t StMCAsymMaker::get_unpolPDF_NLO(int flavor, double x, double Q2){
01420 
01421     Double_t pdf=0.0; 
01422     Int_t iset=1;//NLO MSbar scheme
01423     Int_t er=0;
01424     Int_t fl=10;
01425 
01426     if (flavor==1) fl=2;
01427     if (flavor==2) fl=1;
01428     if (flavor==-1) fl=-2;
01429     if (flavor==-2) fl=-1;
01430     if (flavor==21) fl=0;
01431     if (flavor==3) fl=3;
01432     if (flavor==-3) fl=-3;
01433     if (flavor==4) fl=4;
01434     if (flavor==-4) fl=-4;
01435     if (flavor==5) fl=5;
01436     if (flavor==-5) fl=-5;
01437 
01438     double Q=pow(Q2,0.5);
01439     pdf=ctq5pd_(&iset,&fl,&x,&Q,&er);
01440     if (er!=0) pdf=0.0;
01441 
01442     return pdf;
01443 }
01444 
01445 
01446 Double_t StMCAsymMaker::getPartonicALL(double s, double t, double u, int sub, int inA, int inB, int outA, int outB){
01447 
01448     //Werner definitions:
01449     //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
01450     //PYTHIA definitions:
01451     //1: 11a                       2:11b       3: 12a             4:11 and 12b    5: 13          6: 53         7: 28       8: 68
01452     //NOTES:
01453     // 3==5==6==-1  1==7  1!=2 and 1!=4
01454 
01455     double N1,N2,N3,N4,N5,N6,N7,N8;
01456     double D1,D2,D3,D4,D5,D6,D7,D8;
01457     double all=-10;
01458 
01459     num_(&s,&t,&u,&N1,&N2,&N3,&N4,&N5,&N6,&N7,&N8);
01460     denom_(&s,&t,&u,&D1,&D2,&D3,&D4,&D5,&D6,&D7,&D8);
01461     if (0){
01462         cout<<"s="<<s<<" t="<<t<<" u="<<u<<" sub="<<sub<<" inA="<<inA<<" inB="<<inB<<" outA="<<outA<<" outB="<<outB<<endl;
01463         cout<<" 1="<<N1<<" "<<D1<<endl;
01464         cout<<" 2="<<N2<<" "<<D2<<endl;
01465         cout<<" 3="<<N3<<" "<<D3<<endl;
01466         cout<<" 4="<<N4<<" "<<D4<<endl;
01467         cout<<" 5="<<N5<<" "<<D5<<endl;
01468         cout<<" 6="<<N6<<" "<<D6<<endl;
01469         cout<<" 7="<<N7<<" "<<D7<<endl;
01470         cout<<" 8="<<N8<<" "<<D8<<endl;
01471     }
01472 
01473 
01474     if ((sub==11)&&(abs(inA)!=abs(inB))) all=N1/D1;
01475     if ((sub==11)&&(inA==inB)&&(outA==outB)) all=N2/D2;
01476     //This line added as bug fix 07/17/08 RHF
01477     //Before then these events given all=N2/D2
01478     if ((sub==11)&&(inA==(-1*inB))&&(outA==(-1*outB))&&(inA==outA)&&(inB==outB)) all=N4/D4;
01479     if ((sub==12)&&(abs(inA)!=abs(outA))) all=N3/D3;
01480     if ((sub==12)&&(abs(inA)==abs(outA))) all=N4/D4;
01481     if (sub==13) all=N5/D5;
01482     if (sub==53) all=N6/D6;
01483     if (sub==28) all=N7/D7;
01484     if (sub==68) all=N8/D8;
01485 
01486     //prompt photon subprocesses
01487     if (sub==29) all=N7/D7;//q_i + g -> q_i + gamma
01488     if (sub==14) all=N6/D6;//q_i+qbar_i -> g + gamma
01489     if (sub==18) all=N6/D6;//q_i+qbar_i -> gamma + gamma
01490 
01491     return all;
01492 }
01493 
01494 Double_t StMCAsymMaker::getProtonA1(Double_t x,Double_t Q2){
01495 
01496   // This piece coded by Pibero
01497   // Virtual photon-nucleon asymmetry measured from deep inelatic scattering
01498   // of polarized charged lepton beams from polarized targets. It is defined
01499   // as the sum of polarized PDF weighted by charge square over the sum of
01500   // unpolarized PDF weighted by charge square for all quark/antiquark flavors.
01501 
01502   // Flavors are PDG codes for the quarks and gluons:
01503   // down=1, up=2, strange=3, charm=4, bottom=5, top=6 
01504   // The weights are proportional to the square of the  
01505   // quark charge (in units of e) and indexed by flavor.
01506   // d=-1/3, u=2/3, s=-1/3, c=2/3, b=-1/3, t=2/3.
01507   // The PDG code for an antiquark is just the negative
01508   // of its partner quark. The charge is also the negative
01509   // of its partner quark charge, so the weights are the
01510   // same for quarks and antiquarks. The PDF are GRSV-standard
01511   
01512   const Double_t weights[] = { -1, 1, 4, 1, 4, 1, 4 };
01513   Double_t polSum = 0;
01514   Double_t unpolSum = 0;
01515 
01516   for (int flavor = 1; flavor <= 6; ++flavor) 
01517     {
01518       Double_t polPdf = get_polPDF_NLO(flavor, x, Q2);
01519       polPdf += get_polPDF_NLO(-flavor, x, Q2);
01520       polPdf *= weights[flavor];
01521       polSum += polPdf;
01522       Double_t unpolPdf = get_unpolPDF_NLO(flavor, x, Q2);
01523       unpolPdf += get_unpolPDF_NLO(-flavor, x, Q2);
01524       unpolPdf *= weights[flavor];
01525       unpolSum += unpolPdf;     
01526     }
01527   
01528   Double_t A1 = polSum / unpolSum;
01529   
01530   return A1;
01531 
01532 }
01533 
01534 struct PDFWrapper {
01535   typedef double (*PDF)(int,double,double);
01536 
01537   static const PDF pdfs[StPythiaEvent::NPDF];
01538   static PDF pdf;
01539   static int flavor;
01540   static double Q2;
01541 
01542   static void setPdfFlavorQ2(int pdf, int flavor, double Q2)
01543   {
01544     PDFWrapper::pdf = pdfs[pdf];
01545     PDFWrapper::flavor = flavor;
01546     PDFWrapper::Q2 = Q2;
01547   }
01548 
01549   static double eval(double& x) {  return pdf(flavor,x,Q2); }
01550 };
01551 
01552 const PDFWrapper::PDF PDFWrapper::pdfs[StPythiaEvent::NPDF] = {
01553   StMCAsymMaker::get_polPDF_LO,
01554   StMCAsymMaker::get_polPDF_NLO, // STD
01555   StMCAsymMaker::get_polPDF_NLO_g0,
01556   StMCAsymMaker::get_polPDF_NLO_gmax,
01557   StMCAsymMaker::get_polPDF_NLO_gmin,
01558   StMCAsymMaker::get_polPDF_NLO_m015,
01559   StMCAsymMaker::get_polPDF_NLO_m030,
01560   StMCAsymMaker::get_polPDF_NLO_m045,
01561   StMCAsymMaker::get_polPDF_NLO_m060,
01562   StMCAsymMaker::get_polPDF_NLO_m075,
01563   StMCAsymMaker::get_polPDF_NLO_m090,
01564   StMCAsymMaker::get_polPDF_NLO_m105,
01565   StMCAsymMaker::get_polPDF_NLO_p030,
01566   StMCAsymMaker::get_polPDF_NLO_p045,
01567   StMCAsymMaker::get_polPDF_NLO_p060,
01568   StMCAsymMaker::get_polPDF_NLO_p070,
01569   StMCAsymMaker::get_polPDF_NLO_GSA,
01570   StMCAsymMaker::get_polPDF_NLO_GSB,
01571   StMCAsymMaker::get_polPDF_NLO_GSC,
01572   StMCAsymMaker::get_polPDF_NLO_DSSV,
01573   StMCAsymMaker::get_polPDF_NLO_LSS1,
01574   StMCAsymMaker::get_polPDF_NLO_LSS2,
01575   StMCAsymMaker::get_polPDF_NLO_LSS3,
01576   StMCAsymMaker::get_polPDF_NLO_AAC1,
01577   StMCAsymMaker::get_polPDF_NLO_AAC2,
01578   StMCAsymMaker::get_polPDF_NLO_AAC3,
01579   StMCAsymMaker::get_polPDF_NLO_BB1,
01580   StMCAsymMaker::get_polPDF_NLO_BB2,
01581   StMCAsymMaker::get_polPDF_NLO_DNS1,
01582   StMCAsymMaker::get_polPDF_NLO_DNS2,
01583   StMCAsymMaker::get_polPDF_NLO_DSSV2009a
01584 };
01585 
01586 PDFWrapper::PDF PDFWrapper::pdf = 0;
01587 int PDFWrapper::flavor = 0;
01588 double PDFWrapper::Q2 = 0;
01589 
01590 double StMCAsymMaker::get_polPDF_firstMoment(int pdf, int flavor, double Q2, double xmin, double xmax, double epsilon)
01591 {
01592   PDFWrapper::setPdfFlavorQ2(pdf,flavor,Q2);
01593   return dinteg_(PDFWrapper::eval,xmin,xmax,epsilon);
01594 }