Index: branches/fitting/examples/sartreMinuitMain.cpp =================================================================== --- branches/fitting/examples/sartreMinuitMain.cpp (revision 386) +++ branches/fitting/examples/sartreMinuitMain.cpp (revision 387) @@ -1,2037 +1,1940 @@ #include #include #include #include #include #include #include #include #include #include "AlphaStrong.h" #include "Amplitudes.h" #include "BreakupProduct.h" #include "Constants.h" #include "CrossSection.h" #include "DglapEvolution.h" #include "Enumerations.h" #include "Event.h" #include "Table.h" #include "EventGeneratorSettings.h" #include "Settings.h" #include "TableGeneratorSettings.h" #include "Kinematics.h" #include "InclusiveDiffractiveCrossSections.h" #include "PhotonFlux.h" #include "DipoleModel.h" #include "WaveOverlap.h" #include "DipoleModelParameters.h" #include "Parameters_Minuit2.h" #include "Settings.h" #include "TableGeneratorSettings.h" #include "TLorentzVector.h" #include "Math/WrappedTF1.h" #include "Math/GaussIntegrator.h" #include "TH1D.h" #include "TH2D.h" #include "TFile.h" #include "TF1.h" #include "TGraphErrors.h" #include "TMath.h" #include "Math/Functor.h" #include "Math/IntegratorMultiDim.h" #include "Minuit2/FunctionMinimum.h" #include "Minuit2/MnMigrad.h" #include "Minuit2/MnSimplex.h" #include "Minuit2/MnHesse.h" #include "Minuit2/MnUserParameters.h" #include "Minuit2/MnPrint.h" #include "Minuit2/FCNBase.h" #include "Minuit2/Minuit2Minimizer.h" #include #include #include #include -#include +#include "MinuitGlobalVariables.h" + +#if defined(SARTRE_IN_MULTITHREADED_MODE) +#include +#endif + +// #include #define PR(x) cout << #x << " = " << (x) << endl; using namespace std; double amplitude(const double*); -////////////////////////////////////////////////////////////////////////////////////////////////// -//Global Variables://///////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////////////////////// - -InclusiveDiffractiveCrossSections myIntegrals; -// PhotonFlux myFlux; -//Minuit2Minimizer m2m; - -TableGeneratorSettings* settings = TableGeneratorSettings::instance(); -EventGeneratorSettings* settings2 = EventGeneratorSettings::instance(); -DipoleModelType model = settings->dipoleModelType(); -PhotonFlux mFlux; -DipoleModel_bSat mbSat; -DipoleModel_bNonSat mbNonSat; -WaveOverlapVM mOverlap; -double mVMmass; -double mJPsiMass=3.096916; -double mPhiMass=1.019445; -double mRhoMass=0.775491; -bool mIsTransverse; -double mQ2; -double mT; -double mW2; -double mRelativePrecision=1e-2; -double mS=319*319; -bool mUseRealCorrection=true; -bool mUseSkewednessCorrection=true; -double mRealCorrectionT=1; -double mRealCorrectionL=1; - -DglapEvolution &dglap = DglapEvolution::instance(); - -////////////////////////////////////////////////////////////////////////////////////////////////// -///END GLOBAL VARIABLES/////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////////////////////// - -//////////////////////////////Exclusive Diffraction Plots///////////////////////////////////////// - - - TGraphErrors* gsartreH11 = new TGraphErrors(8); - TGraphErrors* gHERAH11 = new TGraphErrors(8); - TGraphErrors* gsartreH12 = new TGraphErrors(4); - TGraphErrors* gHERAH12 = new TGraphErrors(4); - TGraphErrors* gsartreH13 = new TGraphErrors(4); - TGraphErrors* gHERAH13 = new TGraphErrors(4); - TGraphErrors* gsartreH14 = new TGraphErrors(4); - TGraphErrors* gHERAH14 = new TGraphErrors(4); - - TGraphErrors* gHERAZEUSPhi1 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi2 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi3 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi4 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi5 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi6 = new TGraphErrors(4); - TGraphErrors* gHERAZEUSPhi7 = new TGraphErrors(4); - - TGraphErrors* gsartreZEUSPhi1 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi2 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi3 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi4 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi5 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi6 = new TGraphErrors(4); - TGraphErrors* gsartreZEUSPhi7 = new TGraphErrors(4); - - TGraphErrors* gsartreH11Phi = new TGraphErrors(7); - TGraphErrors* gHERAH11Phi = new TGraphErrors(7); - TGraphErrors* gsartreH12Phi = new TGraphErrors(7); - TGraphErrors* gHERAH12Phi = new TGraphErrors(7); - TGraphErrors* gsartreH13Phi = new TGraphErrors(7); - TGraphErrors* gHERAH13Phi = new TGraphErrors(7); - - TGraphErrors* gsartreZEUS1 = new TGraphErrors(4); - TGraphErrors* gHERAZEUS1 = new TGraphErrors(4); - - TGraphErrors* gsartreZEUS1Phi = new TGraphErrors(7); - TGraphErrors* gHERAZEUS1Phi = new TGraphErrors(7); - - TGraphErrors* gHERAH1JPsi_13HE = new TGraphErrors(6); - TGraphErrors* gHERAH1JPsi_13LE = new TGraphErrors(4); - TGraphErrors* gsartreH1JPsi_13HE = new TGraphErrors(6); - TGraphErrors* gsartreH1JPsi_13LE = new TGraphErrors(4); - - TGraphErrors* gsartreH1Rho1 = new TGraphErrors(7); - TGraphErrors* gsartreH1Rho2 = new TGraphErrors(7); - TGraphErrors* gsartreH1Rho3 = new TGraphErrors(7); - TGraphErrors* gsartreH1Rho4 = new TGraphErrors(7); - TGraphErrors* gsartreH1Rho5 = new TGraphErrors(7); - TGraphErrors* gHERAH1Rho1 = new TGraphErrors(7); - TGraphErrors* gHERAH1Rho2 = new TGraphErrors(7); - TGraphErrors* gHERAH1Rho3 = new TGraphErrors(7); - TGraphErrors* gHERAH1Rho4 = new TGraphErrors(7); - TGraphErrors* gHERAH1Rho5 = new TGraphErrors(7); - -TH1D* hRealCorrectionsT = new TH1D("hRealCorrectionsT", "hRealCorrectionsT", 100, 0., 5.); -TH1D* hRealCorrectionsL = new TH1D("hRealCorrectionsL", "hRealCorrectionsL", 100, 0., 5.); - -////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////// EXCLUSIVE //////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// double dAdx(const double* x){ - double z=x[0]; - double r=x[1]; //fm - double b=x[2]; //fm - double Delta=sqrt(-mT); //GeV - double dsigmadb2p=0; - double dsigmadb2m=0; - - double xp=Kinematics::xprobe(mQ2, mW2, mVMmass); - - double waveOverlap=0; - if(mIsTransverse) waveOverlap=mOverlap.T(z, mQ2, r); //GeV^2 - else waveOverlap=mOverlap.L(z, mQ2, r); //GeV^2 - - double bessel1=TMath::BesselJ0(b*Delta/hbarc); - double bessel2=TMath::BesselJ0((1-z)*r*Delta/hbarc); - - double h=xp*1e-2; - if(model==bSat){ - dsigmadb2p=mbSat.dsigmadb2ep(r, b, xp+h); - dsigmadb2m=mbSat.dsigmadb2ep(r, b, xp-h); - } - if(model==bNonSat){ - dsigmadb2p=mbNonSat.dsigmadb2ep(r, b, xp+h); - dsigmadb2m=mbNonSat.dsigmadb2ep(r, b, xp-h); - } - double derivate=(dsigmadb2p-dsigmadb2m)/(2.*h); - double result=M_PI*r*b/hbarc2*waveOverlap*bessel1*bessel2*derivate; //GeV^0 - - return result; - + double z=x[0]; + double r=x[1]; //fm + double b=x[2]; //fm + double Delta=sqrt(-mT); //GeV + double dsigmadb2p=0; + double dsigmadb2m=0; + + double xp=Kinematics::xprobe(mQ2, mW2, mVMmass); + + double waveOverlap=0; + if(mIsTransverse) waveOverlap=mOverlap.T(z, mQ2, r); //GeV^2 + else waveOverlap=mOverlap.L(z, mQ2, r); //GeV^2 + + double bessel1=TMath::BesselJ0(b*Delta/hbarc); + double bessel2=TMath::BesselJ0((1-z)*r*Delta/hbarc); + + double h=xp*1e-2; + if(model==bSat){ + dsigmadb2p=mbSat.dsigmadb2ep(r, b, xp+h); + dsigmadb2m=mbSat.dsigmadb2ep(r, b, xp-h); + } + if(model==bNonSat){ + dsigmadb2p=mbNonSat.dsigmadb2ep(r, b, xp+h); + dsigmadb2m=mbNonSat.dsigmadb2ep(r, b, xp-h); + } + double derivate=(dsigmadb2p-dsigmadb2m)/(2.*h); + double result=M_PI*r*b/hbarc2*waveOverlap*bessel1*bessel2*derivate; //GeV^0 + + return result; + } double amplitude(const double* x){ - double z=x[0]; - double r=x[1]; //fm - double b=x[2]; //fm - double Delta=sqrt(-mT); //GeV - double dsigmadb2=0; - - double xp=Kinematics::xprobe(mQ2, mW2, mVMmass); - - double waveOverlap=0; - if(mIsTransverse) waveOverlap=mOverlap.T(z, mQ2, r); //GeV^2 - else waveOverlap=mOverlap.L(z, mQ2, r); //GeV^2 - - if(model==bSat) dsigmadb2=mbSat.dsigmadb2ep(r, b, xp); - if(model==bNonSat) dsigmadb2=mbNonSat.dsigmadb2ep(r, b, xp); - - double bessel1=TMath::BesselJ0(b*Delta/hbarc); - double bessel2=TMath::BesselJ0((1-z)*r*Delta/hbarc); - - double result=M_PI*r*b/hbarc2*waveOverlap*bessel1*bessel2*dsigmadb2; //GeV^0 - - return result; + double z=x[0]; + double r=x[1]; //fm + double b=x[2]; //fm + double Delta=sqrt(-mT); //GeV + double dsigmadb2=0; + + double xp=Kinematics::xprobe(mQ2, mW2, mVMmass); + + double waveOverlap=0; + if(mIsTransverse) waveOverlap=mOverlap.T(z, mQ2, r); //GeV^2 + else waveOverlap=mOverlap.L(z, mQ2, r); //GeV^2 + + if(model==bSat) dsigmadb2=mbSat.dsigmadb2ep(r, b, xp); + if(model==bNonSat) dsigmadb2=mbNonSat.dsigmadb2ep(r, b, xp); + + double bessel1=TMath::BesselJ0(b*Delta/hbarc); + double bessel2=TMath::BesselJ0((1-z)*r*Delta/hbarc); + + double result=M_PI*r*b/hbarc2*waveOverlap*bessel1*bessel2*dsigmadb2; //GeV^0 + + return result; } double exclusiveCS(const double* x){ - mQ2=x[0]; - mW2=x[1]; - mT=-x[2]; - - double rbmax=sqrt(2)*hbarc*6; - double a[3]={0., 0., 0.}; //z ,r, b - double b[3]={1., rbmax, rbmax}; - ROOT::Math::Functor amp(&litude, 3); - ROOT::Math::IntegratorMultiDim igAmp(ROOT::Math::IntegrationMultiDim::kADAPTIVE); - igAmp.SetFunction(amp); - igAmp.SetAbsTolerance(0.); - igAmp.SetRelTolerance(mRelativePrecision/10.); - - double amps[2]={0., 0.}; - for(int i=0; i<2; i++){ - mIsTransverse=i; - - amps[i]=igAmp.Integral(a, b); //fm2 - mRealCorrectionT=1.; - mRealCorrectionL=1.; - if(mUseRealCorrection){ - double h=mW2*1e-5; - mW2+=h; - double amp_p=igAmp.Integral(a, b); - mW2-=2*h; - double amp_m=igAmp.Integral(a, b); - mW2+=h; - - - double dlogAdW2=log(amp_p/amp_m)/(2.*h); - //dW2/dlog(1/x): - double dW2dlogOneOverxp=mW2-protonMass2+mQ2; - //lambda=dlog(A)/dlog(1/x)=dlog(A)/dW2*dW2/dlog(1/x): - double lambda = dlogAdW2*dW2dlogOneOverxp; - double beta=tan(M_PI/2.*lambda); - if(beta*beta>0.5) beta=1./sqrt(2); - if(mIsTransverse){ - mRealCorrectionT=(1+beta*beta); - hRealCorrectionsT->Fill(mRealCorrectionT); - } - else { - mRealCorrectionL=(1+beta*beta); - hRealCorrectionsL->Fill(mRealCorrectionL); - } - if(1+beta*beta>2.) cout<<"corr="<<1+beta*beta<<", lambda="<> qsqr; l >> x; l >> y; l >> sigmar; l >> err; - - if (StrToReal(x)>maxx or StrToReal(x)maxQ2) continue; - points++; - - //MZ x cuts - double xbj=StrToReal(x)*(1.+4.*pow(params_xtil_c_mass, 2)/StrToReal(qsqr)); - if(xbj>=1e-2){ - continue; - } - - Q2vals.push_back(StrToReal(qsqr)); - xvals.push_back(StrToReal(x)); - yvals.push_back(StrToReal(y)); - sigmarvals.push_back(StrToReal(sigmar)); - errors.push_back(StrToReal(err)); - only_charm.push_back(onlycharm); - } - - cout << "# Loaded " << points << " datapoints from " << filename << " in Q2 range " << minQ2 << " - " << maxQ2 << " GeV^2, now we have in total " << sigmarvals.size() << " points " << endl; - - return 0; + + ifstream file(filename.c_str()); + if(file.fail()){ + cout<<"File does not exist:"<> qsqr; l >> x; l >> y; l >> sigmar; l >> err; + + if (StrToReal(x)>maxx or StrToReal(x)maxQ2) continue; + points++; + + //MZ x cuts + double xbj=StrToReal(x)*(1.+4.*pow(params_xtil_c_mass, 2)/StrToReal(qsqr)); + if(xbj>=1e-2){ + continue; + } + // cout<<"Q2="<= NumPoints()) - { - cerr << "Index " << n << " out of range! " << endl; - exit(1); - } - return Q2vals[n]; + if (n>= NumPoints()) + { + cerr << "Index " << n << " out of range! " << endl; + exit(1); + } + return Q2vals[n]; } double Data::xVals(unsigned int n) const { - return xvals[n]; + return xvals[n]; } double Data::yVals(unsigned int n) const { - return yvals[n]; + return yvals[n]; } double Data::sigmar(unsigned int n) const { - return sigmarvals[n]; + return sigmarvals[n]; } double Data::sigmaErr(unsigned int n) const { - return errors[n]; + return errors[n]; } bool Data::OnlyCharm(unsigned int n) const { - return only_charm[n]; + return only_charm[n]; } Data::Data() { - minx = 1e-99; - maxx = 1e-2; - minQ2 = 1.5; //BS MZ sets minQ2 at 1.5; // was 0.001; - maxQ2 = 50; //BS MZ sets maxQ2 at 50; // was 650; - weight = 1.0; //no weights at the moment + minx = 1e-99; + maxx = 1e-2; + minQ2 = 0.;//1.5; //BS MZ sets minQ2 at 1.5; // was 0.001; + maxQ2 = 50.;//50; //BS MZ sets maxQ2 at 50; // was 650; + weight = 1.0; //no weights at the moment } double uiSigma(const double* var){ - double Q2=var[0]; - double W2=var[1]; - double t=var[2]; - double z=var[3]; - double b=var[4]; - double r=var[5]; - double xp=Kinematics::xprobe(Q2, W2, mVMmass); - double Delta=sqrt(-t); - - double waveOverlap=mOverlap.T(z, Q2, r); - double dipole=mbSat.dsigmadb2ep(r, b, xp); - double BesselJ0r=TMath::BesselJ0((1-z)*r/hbarc*Delta); - double BesselJ0b=TMath::BesselJ0(b/hbarc*Delta); - - double amplitude=M_PI*r/hbarc*waveOverlap*dipole*BesselJ0r*BesselJ0b*b/hbarc; - - double sigma=1./16./M_PI*amplitude*amplitude; - // PR(Q2); - // PR(W2); - // PR(t); - // PR(z); - // PR(b); - // PR(r); - // PR(xp); - // PR(Delta); - // PR(BesselJ0r); - // PR(BesselJ0b); - // PR(waveOverlap); - // PR(dipole); - // PR(amplitude); - // PR(sigma); - - return sigma; - + double Q2=var[0]; + double W2=var[1]; + double t=var[2]; + double z=var[3]; + double b=var[4]; + double r=var[5]; + double xp=Kinematics::xprobe(Q2, W2, mVMmass); + double Delta=sqrt(-t); + + double waveOverlap=mOverlap.T(z, Q2, r); + double dipole=mbSat.dsigmadb2ep(r, b, xp); + double BesselJ0r=TMath::BesselJ0((1-z)*r/hbarc*Delta); + double BesselJ0b=TMath::BesselJ0(b/hbarc*Delta); + + double amplitude=M_PI*r/hbarc*waveOverlap*dipole*BesselJ0r*BesselJ0b*b/hbarc; + + double sigma=1./16./M_PI*amplitude*amplitude; + // PR(Q2); + // PR(W2); + // PR(t); + // PR(z); + // PR(b); + // PR(r); + // PR(xp); + // PR(Delta); + // PR(BesselJ0r); + // PR(BesselJ0b); + // PR(waveOverlap); + // PR(dipole); + // PR(amplitude); + // PR(sigma); + + return sigma; + } double uiRAve(const double* var){ - double sigma=uiSigma(var); - double r=var[5]; - //PR(sigma); - return r/hbarc*sigma; - + double sigma=uiSigma(var); + double r=var[5]; + //PR(sigma); + return r/hbarc*sigma; + } double calculateMeanR(){ - double down[6]={0.05, 40*40, -1., 0., 0., 0.};// Q2, W2, t, z, b, r - double up[6]={22.4, 160*160, 0., 1., 10., 10.}; - - ROOT::Math::Functor wfs(&uiSigma, 6); - ROOT::Math::IntegratorMultiDim igs(ROOT::Math::IntegrationMultiDim::kADAPTIVE); - igs.SetFunction(wfs); - igs.SetAbsTolerance(0.); - igs.SetRelTolerance(1e-3); - double sigma=igs.Integral(down, up)/hbarc2; - //PR(sigma); - - ROOT::Math::Functor wfr(&uiRAve, 6); - ROOT::Math::IntegratorMultiDim igr(ROOT::Math::IntegrationMultiDim::kADAPTIVE); - igr.SetFunction(wfr); - igr.SetAbsTolerance(0.); - igr.SetRelTolerance(1e-3); - double rsigma=igr.Integral(down, up)/hbarc2; - //PR(rsigma); - - double result=rsigma/sigma; //GeV-1 - //PR(result); - - return result; + double down[6]={0.05, 40*40, -1., 0., 0., 0.};// Q2, W2, t, z, b, r + double up[6]={22.4, 160*160, 0., 1., 10., 10.}; + + ROOT::Math::Functor wfs(&uiSigma, 6); + ROOT::Math::IntegratorMultiDim igs(ROOT::Math::IntegrationMultiDim::kADAPTIVE); + igs.SetFunction(wfs); + igs.SetAbsTolerance(0.); + igs.SetRelTolerance(1e-3); + double sigma=igs.Integral(down, up)/hbarc2; + //PR(sigma); + + ROOT::Math::Functor wfr(&uiRAve, 6); + ROOT::Math::IntegratorMultiDim igr(ROOT::Math::IntegrationMultiDim::kADAPTIVE); + igr.SetFunction(wfr); + igr.SetAbsTolerance(0.); + igr.SetRelTolerance(1e-3); + double rsigma=igr.Integral(down, up)/hbarc2; + //PR(rsigma); + + double result=rsigma/sigma; //GeV-1 + //PR(result); + + return result; } void DISFitter::AddDataset(Data& data){ - datasets.push_back(&data); + datasets.push_back(&data); } DISFitter::DISFitter(MnUserParameters params){ - parameters = params; + parameters = params; } double DISFitter::operator()(const std::vector& par) const { - // arg par - // 0 BG - // 1 m_uds - // 2 m_C - // 3 m_B - // 4 C - // 5 mu02 - // 6 Ag - // 7 lambdaG - // 8 Rshrink - // 9 as_qqg - // 10 Rhs (BG) - // 11 nLP (BG) - // 12 R0WS (n/a) - // 13 dWS (n/a) - // 14 uCE (BG) - // 15 EP (??) - // 16 BV (??) - // 17 RV (??) - // 18 wV (??) - - - cout << "----------- Iteration: " << runCount << " ------------" << endl; - runCount++; - - if(par[5]!=temp[5] || par[6]!=temp[6] || par[7]!=temp[7]) - { - compDGLAP=true; - cout<<"DGLAP changed"< parametersForSettings; - parametersForSettings.resize(12); - parametersForSettings={par[1], par[1], par[1], par[2], par[3], //quarkMasses - par[0], //BG - par[5], //mu02 - par[7], //lambdaG - par[6], //Ag - par[4], //C - params_Shape, - 2}; //CE V - settings->setDipoleModelCustomParameters(parametersForSettings); - - FitParameters fitparams; - fitparams.values = ∥ - fitparams.parameter = ¶meters; - - - if(compDGLAP==true){ - // - // DGLAP Evolution can be speed up by using lookup tables - // - dglap.reinit(); - dglap.generateLookupTable(1000, 1000); - dglap.useLookupTable(true); - } - else { - cout<< "DGLAP parameters unchanged"<NumPoints(); - - //int maxPoints=data.NumPoints(); - double chi=0; - double chi2=0; - int ndf=0; - - //Data inclusion switches - bool HERA_I_DIS=true; - bool Charm=true; - bool combinedQ2overbeta=false; - bool H1_JPsi=false; - bool ZEUS_JPsi_electro=H1_JPsi; - bool H1_JPsi_13=H1_JPsi; - bool H1_Rho=false; - bool H1_Phi=false; - bool ZEUS_Phi_electro=H1_Phi; - - - if(HERA_I_DIS){ - PR(HERA_I_DIS); - for(unsigned int dataset=0; datasetNumPoints(); i++){ - //for(int i = 0; i<10; i++){ - //double sqrts=datasets[dataset]->sqrtsVals(i); - //cout<<"sqrts= "<Q2Vals(i);//data.Q2vals[i]; - double x=datasets[dataset]->xVals(i);//data.xvals[i]; - double y=datasets[dataset]->yVals(i);//data.yvals[i]; //replace with Kinematics::y for diffraction data - //double xbj=x*(1.+4.*pow(params_xtil_c_mass, 2)/Q2); - // if(xbj>=1e-2){ - // cout<<"xtil out of range"<yVals(i);//data.yvals[i]; //replace with Kinematics::y for diffraction data - double W2=Kinematics::W2(Q2, x); - - //cout<<"DIS"<sigmar(i))/(datasets[dataset]->sigmaErr(i)); - //double chipnt=(sigma_r-data.sigmarvals[i])/data.errors[i]; - // #pragma omp critical - // { - - //for SK: add print statements at chi, chi2, kinematics, calculateTotal -> THREADWISE - //using omp_num_threads variable - double data=datasets[dataset]->sigmar(i); - double chipnt=(sigma_r-data)/(datasets[dataset]->sigmaErr(i)); - chi+=chipnt; - //PR(chipnt); - chi2+=datasets[dataset]->Weight()*chipnt*chipnt; - // PR(chipnt); - // PR(sigma_r); - // PR(data); - // PR(sigma_r/data); - //chi2+= ((sigma_r-data.sigmarvals[i])/data.errors[i])*((sigma_r-data.sigmarvals[i])/data.errors[i]); - ndf+=datasets[dataset]->Weight(); - - - } - //cout << "Dataset: " << dataset << "loop count: " << i << "----" << endl; - } - cout<<"======================================================================="<sqrtsVals(i); - //cout<<"sqrts= "<maxQ2 || Q2maxx || x=1e-2){ - //cout<<"xtil out of range"<yVals(i);//data.yvals[i]; //replace with Kinematics::y for diffraction data - double W2=Kinematics::W2(Q2, x); - - // #pragma omp critical - { - myIntegrals.setWhichFlavours(false, false, false, true, false); - - //cout<<"CHARM"<sigmar(i))/(datasets[dataset]->sigmaErr(i)); - //double chipnt=(sigma_r-data.sigmarvals[i])/data.errors[i]; - // #pragma omp critical - // { - - //for SK: add print statements at chi, chi2, kinematics, calculateTotal -> THREADWISE - //using omp_num_threads variable - - chi+=(sigma_r-sigmar[i])/(sigmaErr[i]); - chi2+=weight*((sigma_r-sigmar[i])/(sigmaErr[i]))* - ((sigma_r-sigmar[i])/(sigmaErr[i])); - //chi2+= ((sigma_r-data.sigmarvals[i])/data.errors[i])*((sigma_r-data.sigmarvals[i])/data.errors[i]); - ndf+=weight; - } - //cout << "Dataset: " << dataset << "loop count: " << i << "----" << endl; - - cout<<"======================================================================="<CORR is "<nb - - TGraphErrors *gsartre1 = new TGraphErrors(4); - TGraphErrors *gHERAData1 = new TGraphErrors(4); - - for(int i=0; i<4; i++){ - double beta=beta1[i]; - if(betabetaMax) continue; - double Q2=Q21[i]; - double x=xpom*beta; - double W2=Kinematics::W2(Q2, x); - myIntegrals.calculate_dsigmadbeta_QQ(Q2, W2, beta); - if(betaMin<0.5){ - myIntegrals.calculate_dsigmadbeta_QQG(Q2, W2, beta); //BS QQG - } - double fact=Q2*beta/(4*M_PI*M_PI*alpha_em); //dsigma/dbeta-->xF - double xFT=0; - double xFL=0; - for(int j=0; jSetPoint(i, variable, xsigma_r); - gHERAData1->SetPoint(i, variable, HERA1[i]); - gHERAData1->SetPointError(i, 0, err1[i]); - - double chipnt=(xsigma_r-HERA1[i])/err1[i]; - PR(chipnt); - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - } - PR(chi); - cout<<"chi2="<Write("gsartre1"); - // gHERAData1->Write("gHERAData1"); - // gsartre2->Write("gsartre2"); - // gHERAData2->Write("gHERAData2"); - // gsartre3->Write("gsartre3"); - // gHERAData3->Write("gHERAData3"); - // gsartre4->Write("gsartre4"); - // gHERAData4->Write("gHERAData4"); - - // gsartre5->Write("gsartre5"); - // gHERAData5->Write("gHERAData5"); - // hfile->Close(); - // cout << rootfile.c_str() <<" written." << endl; - - cout<<"======================================================================="<SetPoint(i, tvals1[i], H11[i]); - gHERAH11->SetPointError(i, 0., error1[i]); - gsartreH11->SetPoint(i, tvals1[i], sigma); - - double chipnt=(sigma-H11[i])/error1[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals2[i], H12[i]); - gHERAH12->SetPointError(i, 0., error2[i]); - gsartreH12->SetPoint(i, tvals2[i], sigma); - - double chipnt=(sigma-H12[i])/error2[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals2[i], H13[i]); - gHERAH13->SetPointError(i, 0., error3[i]); - gsartreH13->SetPoint(i, tvals2[i], sigma); - - double chipnt=(sigma-H13[i])/error3[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals2[i], H14[i]); - gHERAH14->SetPointError(i, 0., error4[i]); - gsartreH14->SetPoint(i, tvals2[i], sigma); - - double chipnt=(sigma-H14[i])/error4[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], ZEUS1[i]); - gHERAZEUS1->SetPointError(i, 0., error1[i]); - gsartreZEUS1->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-ZEUS1[i])/error1[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["< parametersForSettings; + parametersForSettings.resize(12); + parametersForSettings={par[1], par[1], par[1], par[2], par[3], //quarkMasses + par[0], //BG + par[5], //mu02 + par[7], //lambdaG + par[6], //Ag + par[4], //C + params_Shape, + 2}; //CE V + settings->setDipoleModelCustomParameters(parametersForSettings); + + FitParameters fitparams; + fitparams.values = ∥ + fitparams.parameter = ¶meters; + + + if(compDGLAP==true){ + // + // DGLAP Evolution can be speed up by using lookup tables + // + dglap.reinit(); + dglap.generateLookupTable(1000, 1000); + dglap.useLookupTable(true); + } + else { + cout<< "DGLAP parameters unchanged"<NumPoints(); + + //int maxPoints=data.NumPoints(); + double chi=0; + double chi2=0; + int ndf=0; + + //Data inclusion switches + bool HERA_I_DIS=true; + bool Charm=true; + bool combinedQ2overbeta=false; + bool H1_JPsi=false; + bool ZEUS_JPsi_electro=H1_JPsi; + bool H1_JPsi_13=H1_JPsi; + bool H1_Rho=false; + bool H1_Phi=false; + bool ZEUS_Phi_electro=H1_Phi; + + + myIntegrals.setParamsDCS(par); + myIntegrals.dipoleModel()->createNormTable(); + if(HERA_I_DIS){ +// int numCPU = boost::thread::hardware_concurrency(); +//#if defined(SARTRE_IN_MULTITHREADED_MODE) // multithreaded +// vector vIntegrals; + // + // Create a vector containing the threads: + // +// std::vector vThreads; +// vThreads.clear(); + // + // Create the thread group: + // +// boost::thread_group gThreads; +// vThreads.push_back(new boost::thread(boost::ref(*mIntegrals.at(i)), +// t, Q2, W2)); +// +// gThreads.add_thread(vThreads.at(i)); + + +//#else // unforked version + for(unsigned int dataset=0; datasetNumPoints(); i++){ + double Q2=datasets[dataset]->Q2Vals(i); + double x=datasets[dataset]->xVals(i); + double y=datasets[dataset]->yVals(i); + double W2=Kinematics::W2(Q2, x); + if(dataset==0) { + myIntegrals.setWhichFlavours(true, true, true, true, true); + } + + myIntegrals.operator()(Q2, W2); + double sigmaT=myIntegrals.mTotalCST; + double sigmaL=myIntegrals.mTotalCSL; + double factor=Q2/(4*M_PI*M_PI*alpha_em); + double F2=factor*(sigmaT+sigmaL); + double FL=factor*sigmaL; + double sigma_r=F2-y*y/(1+(1-y)*(1-y))*FL; + double data=datasets[dataset]->sigmar(i); + double chipnt=(sigma_r-data)/(datasets[dataset]->sigmaErr(i)); + sartreSum+=sigma_r; + dataSum+=data; + chi+=chipnt; + chi2+=datasets[dataset]->Weight()*chipnt*chipnt; + ndf+=datasets[dataset]->Weight(); +// cout<<"sigma_r("<sigmaErr(i)<<", chipnt="<sqrtsVals(i); + //cout<<"sqrts= "<maxQ2 || Q2maxx || x=1e-2){ + //cout<<"xtil out of range"<yVals(i);//data.yvals[i]; //replace with Kinematics::y for diffraction data + double W2=Kinematics::W2(Q2, x); + + // #pragma omp critical + { + myIntegrals.setWhichFlavours(false, false, false, true, false); + + //cout<<"CHARM"<sigmar(i))/(datasets[dataset]->sigmaErr(i)); + //double chipnt=(sigma_r-data.sigmarvals[i])/data.errors[i]; + // #pragma omp critical + // { + + //for SK: add print statements at chi, chi2, kinematics, calculateTotal -> THREADWISE + //using omp_num_threads variable + + chi+=(sigma_r-sigmar[i])/(sigmaErr[i]); + chi2+=weight*((sigma_r-sigmar[i])/(sigmaErr[i]))* + ((sigma_r-sigmar[i])/(sigmaErr[i])); + //chi2+= ((sigma_r-data.sigmarvals[i])/data.errors[i])*((sigma_r-data.sigmarvals[i])/data.errors[i]); + ndf+=weight; + } + //cout << "Dataset: " << dataset << "loop count: " << i << "----" << endl; + + cout<<"======================================================================="<CORR is "<nb + + TGraphErrors *gsartre1 = new TGraphErrors(4); + TGraphErrors *gHERAData1 = new TGraphErrors(4); + + for(int i=0; i<4; i++){ + double beta=beta1[i]; + if(betabetaMax) continue; + double Q2=Q21[i]; + double x=xpom*beta; + double W2=Kinematics::W2(Q2, x); + myIntegrals.calculate_dsigmadbeta_QQ(Q2, W2, beta); + if(betaMin<0.5){ + myIntegrals.calculate_dsigmadbeta_QQG(Q2, W2, beta); //BS QQG + } + double fact=Q2*beta/(4*M_PI*M_PI*alpha_em); //dsigma/dbeta-->xF + double xFT=0; + double xFL=0; + for(int j=0; jSetPoint(i, variable, xsigma_r); + gHERAData1->SetPoint(i, variable, HERA1[i]); + gHERAData1->SetPointError(i, 0, err1[i]); + + double chipnt=(xsigma_r-HERA1[i])/err1[i]; + PR(chipnt); + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + } + PR(chi); + cout<<"chi2="<Write("gsartre1"); + // gHERAData1->Write("gHERAData1"); + // gsartre2->Write("gsartre2"); + // gHERAData2->Write("gHERAData2"); + // gsartre3->Write("gsartre3"); + // gHERAData3->Write("gHERAData3"); + // gsartre4->Write("gsartre4"); + // gHERAData4->Write("gHERAData4"); + + // gsartre5->Write("gsartre5"); + // gHERAData5->Write("gHERAData5"); + // hfile->Close(); + // cout << rootfile.c_str() <<" written." << endl; + + cout<<"======================================================================="<SetPoint(i, tvals1[i], H11[i]); + gHERAH11->SetPointError(i, 0., error1[i]); + gsartreH11->SetPoint(i, tvals1[i], sigma); + + double chipnt=(sigma-H11[i])/error1[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals2[i], H12[i]); + gHERAH12->SetPointError(i, 0., error2[i]); + gsartreH12->SetPoint(i, tvals2[i], sigma); + + double chipnt=(sigma-H12[i])/error2[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals2[i], H13[i]); + gHERAH13->SetPointError(i, 0., error3[i]); + gsartreH13->SetPoint(i, tvals2[i], sigma); + + double chipnt=(sigma-H13[i])/error3[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals2[i], H14[i]); + gHERAH14->SetPointError(i, 0., error4[i]); + gsartreH14->SetPoint(i, tvals2[i], sigma); + + double chipnt=(sigma-H14[i])/error4[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], ZEUS1[i]); + gHERAZEUS1->SetPointError(i, 0., error1[i]); + gsartreZEUS1->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-ZEUS1[i])/error1[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvalsHE[i], sigma); - gHERAH1JPsi_13HE->SetPoint(i, tvalsHE[i], H1valsHE[i]); - gHERAH1JPsi_13HE->SetPointError(i, 0, H1errsHE[i]); - - double chipnt=(sigma-H1valsHE[i])/H1errsHE[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvalsLE[i], sigma); - gHERAH1JPsi_13LE->SetPoint(i, tvalsLE[i], H1valsLE[i]); - gHERAH1JPsi_13LE->SetPointError(i, 0, H1errsLE[i]); - - double chipnt=(sigma-H1valsLE[i])/H1errsLE[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvalsHE[i], sigma); + gHERAH1JPsi_13HE->SetPoint(i, tvalsHE[i], H1valsHE[i]); + gHERAH1JPsi_13HE->SetPointError(i, 0, H1errsHE[i]); + + double chipnt=(sigma-H1valsHE[i])/H1errsHE[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvalsLE[i], sigma); + gHERAH1JPsi_13LE->SetPoint(i, tvalsLE[i], H1valsLE[i]); + gHERAH1JPsi_13LE->SetPointError(i, 0, H1errsLE[i]); + + double chipnt=(sigma-H1valsLE[i])/H1errsLE[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], H1vals1[i]); - gHERAH1Rho1->SetPointError(i, 0., H1errs1[i]); - gsartreH1Rho1->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-H1vals1[i])/H1errs1[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], H1vals2[i]); - gHERAH1Rho2->SetPointError(i, 0., H1errs2[i]); - gsartreH1Rho2->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-H1vals2[i])/H1errs2[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], H1vals3[i]); - gHERAH1Rho3->SetPointError(i, 0., H1errs3[i]); - gsartreH1Rho3->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-H1vals3[i])/H1errs3[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], H1vals4[i]); - gHERAH1Rho4->SetPointError(i, 0., H1errs4[i]); - gsartreH1Rho4->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-H1vals4[i])/H1errs4[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], H1vals5[i]); - gHERAH1Rho5->SetPointError(i, 0., H1errs5[i]); - gsartreH1Rho5->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-H1vals5[i])/H1errs5[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], H1vals1[i]); + gHERAH1Rho1->SetPointError(i, 0., H1errs1[i]); + gsartreH1Rho1->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-H1vals1[i])/H1errs1[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], H1vals2[i]); + gHERAH1Rho2->SetPointError(i, 0., H1errs2[i]); + gsartreH1Rho2->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-H1vals2[i])/H1errs2[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], H1vals3[i]); + gHERAH1Rho3->SetPointError(i, 0., H1errs3[i]); + gsartreH1Rho3->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-H1vals3[i])/H1errs3[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], H1vals4[i]); + gHERAH1Rho4->SetPointError(i, 0., H1errs4[i]); + gsartreH1Rho4->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-H1vals4[i])/H1errs4[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], H1vals5[i]); + gHERAH1Rho5->SetPointError(i, 0., H1errs5[i]); + gsartreH1Rho5->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-H1vals5[i])/H1errs5[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, t, sigma); - gHERAH11Phi->SetPoint(i, t, H11[i]); - gHERAH11Phi->SetPointError(i, 0, error1[i]); - - double chipnt=(sigma-H11[i])/error1[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, t, sigma); + gHERAH11Phi->SetPoint(i, t, H11[i]); + gHERAH11Phi->SetPointError(i, 0, error1[i]); + + double chipnt=(sigma-H11[i])/error1[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], ZEUS1[i]); - gHERAZEUSPhi1->SetPointError(i, 0., error1[i]); - gsartreZEUSPhi1->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-ZEUS1[i])/error1[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], ZEUS2[i]); - gHERAZEUSPhi2->SetPointError(i, 0., error2[i]); - gsartreZEUSPhi2->SetPoint(i, tvals[i], sigma); - - double chipnt=(sigma-ZEUS2[i])/error2[i]; - chi+=chipnt; - chi2+=chipnt*chipnt; - ndf++; - cout<<"chi["<SetPoint(i, tvals[i], ZEUS3[i]); - // gHERAZEUSPhi3->SetPointError(i, 0., error3[i]); - // gsartreZEUSPhi3->SetPoint(i, tvals[i], sigma); - - // double chipnt=(sigma-ZEUS3[i])/error3[i]; - // chi+=chipnt; - // chi2+=chipnt*chipnt; - // ndf++; - // cout<<"chi["<SetPoint(i, tvals[i], ZEUS4[i]); - // gHERAZEUSPhi4->SetPointError(i, 0., error4[i]); - // gsartreZEUSPhi4->SetPoint(i, tvals[i], sigma); - - // double chipnt=(sigma-ZEUS4[i])/error4[i]; - // chi+=chipnt; - // chi2+=chipnt*chipnt; - // ndf++; - // cout<<"chi["<SetPoint(i, tvals[i], ZEUS5[i]); - // gHERAZEUSPhi5->SetPointError(i, 0., error5[i]); - // gsartreZEUSPhi5->SetPoint(i, tvals[i], sigma); - - // double chipnt=(sigma-ZEUS5[i])/error5[i]; - // chi+=chipnt; - // chi2+=chipnt*chipnt; - // ndf++; - // cout<<"chi["<SetPoint(i, tvals[i], ZEUS6[i]); - // gHERAZEUSPhi6->SetPointError(i, 0., error6[i]); - // gsartreZEUSPhi6->SetPoint(i, tvals[i], sigma); - - // double chipnt=(sigma-ZEUS6[i])/error6[i]; - // chi+=chipnt; - // chi2+=chipnt*chipnt; - // ndf++; - // cout<<"chi["<SetPoint(i, tvals[i], ZEUS7[i]); - // gHERAZEUSPhi7->SetPointError(i, 0., error7[i]); - // gsartreZEUSPhi7->SetPoint(i, tvals[i], sigma); - - // double chipnt=(sigma-ZEUS7[i])/error7[i]; - // chi+=chipnt; - // chi2+=chipnt*chipnt; - // ndf++; - // cout<<"chi["<SetPoint(i, tvals[i], ZEUS1[i]); + gHERAZEUSPhi1->SetPointError(i, 0., error1[i]); + gsartreZEUSPhi1->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-ZEUS1[i])/error1[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], ZEUS2[i]); + gHERAZEUSPhi2->SetPointError(i, 0., error2[i]); + gsartreZEUSPhi2->SetPoint(i, tvals[i], sigma); + + double chipnt=(sigma-ZEUS2[i])/error2[i]; + chi+=chipnt; + chi2+=chipnt*chipnt; + ndf++; + cout<<"chi["<SetPoint(i, tvals[i], ZEUS3[i]); + // gHERAZEUSPhi3->SetPointError(i, 0., error3[i]); + // gsartreZEUSPhi3->SetPoint(i, tvals[i], sigma); + + // double chipnt=(sigma-ZEUS3[i])/error3[i]; + // chi+=chipnt; + // chi2+=chipnt*chipnt; + // ndf++; + // cout<<"chi["<SetPoint(i, tvals[i], ZEUS4[i]); + // gHERAZEUSPhi4->SetPointError(i, 0., error4[i]); + // gsartreZEUSPhi4->SetPoint(i, tvals[i], sigma); + + // double chipnt=(sigma-ZEUS4[i])/error4[i]; + // chi+=chipnt; + // chi2+=chipnt*chipnt; + // ndf++; + // cout<<"chi["<SetPoint(i, tvals[i], ZEUS5[i]); + // gHERAZEUSPhi5->SetPointError(i, 0., error5[i]); + // gsartreZEUSPhi5->SetPoint(i, tvals[i], sigma); + + // double chipnt=(sigma-ZEUS5[i])/error5[i]; + // chi+=chipnt; + // chi2+=chipnt*chipnt; + // ndf++; + // cout<<"chi["<SetPoint(i, tvals[i], ZEUS6[i]); + // gHERAZEUSPhi6->SetPointError(i, 0., error6[i]); + // gsartreZEUSPhi6->SetPoint(i, tvals[i], sigma); + + // double chipnt=(sigma-ZEUS6[i])/error6[i]; + // chi+=chipnt; + // chi2+=chipnt*chipnt; + // ndf++; + // cout<<"chi["<SetPoint(i, tvals[i], ZEUS7[i]); + // gHERAZEUSPhi7->SetPointError(i, 0., error7[i]); + // gsartreZEUSPhi7->SetPoint(i, tvals[i], sigma); + + // double chipnt=(sigma-ZEUS7[i])/error7[i]; + // chi+=chipnt; + // chi2+=chipnt*chipnt; + // ndf++; + // cout<<"chi["<rootfile(); TFile *hfile = 0; hfile = new TFile(rootfile.c_str(), "RECREATE"); - + gsartreH11->Write("gsartreH11"); gHERAH11->Write("gHERAH11"); gsartreH12->Write("gsartreH12"); gHERAH12->Write("gHERAH12"); gsartreH13->Write("gsartreH13"); gHERAH13->Write("gHERAH13"); gsartreH14->Write("gsartreH14"); gHERAH14->Write("gHERAH14"); gsartreH11Phi->Write("gsartreH11Phi"); gHERAH11Phi->Write("gHERAH11Phi"); gsartreH12Phi->Write("gsartreH12Phi"); gHERAH12Phi->Write("gHERAH12Phi"); gsartreH13Phi->Write("gsartreH13Phi"); gHERAH13Phi->Write("gHERAH13Phi"); gsartreZEUS1->Write("gsartreZEUS1"); gHERAZEUS1->Write("gHERAZEUS1"); gsartreZEUS1->Write("gsartreZEUS1Phi"); gHERAZEUS1->Write("gHERAZEUS1Phi"); - + gHERAZEUSPhi1->Write("gHERAZEUSPhi1"); gHERAZEUSPhi2->Write("gHERAZEUSPhi2"); gHERAZEUSPhi3->Write("gHERAZEUSPhi3"); gHERAZEUSPhi4->Write("gHERAZEUSPhi4"); gHERAZEUSPhi5->Write("gHERAZEUSPhi5"); gHERAZEUSPhi6->Write("gHERAZEUSPhi6"); gHERAZEUSPhi7->Write("gHERAZEUSPhi7"); - + gsartreZEUSPhi1->Write("gsartreZEUSPhi1"); gsartreZEUSPhi2->Write("gsartreZEUSPhi2"); gsartreZEUSPhi3->Write("gsartreZEUSPhi3"); gsartreZEUSPhi4->Write("gsartreZEUSPhi4"); gsartreZEUSPhi5->Write("gsartreZEUSPhi5"); gsartreZEUSPhi6->Write("gsartreZEUSPhi6"); gsartreZEUSPhi7->Write("gsartreZEUSPhi7"); - + gsartreH1JPsi_13HE->Write("gsartreH1JPsi_13HE"); gHERAH1JPsi_13HE->Write("gHERAH1JPsi_13HE"); gsartreH1JPsi_13LE->Write("gsartreH1JPsi_13LE"); gHERAH1JPsi_13LE->Write("gHERAH1JPsi_13LE"); - + gsartreH1Rho1->Write("gsartreH1Rho1"); gsartreH1Rho2->Write("gsartreH1Rho2"); gsartreH1Rho3->Write("gsartreH1Rho3"); gsartreH1Rho4->Write("gsartreH1Rho4"); gsartreH1Rho5->Write("gsartreH1Rho5"); - + gHERAH1Rho1->Write("gHERAH1Rho1"); gHERAH1Rho2->Write("gHERAH1Rho2"); gHERAH1Rho3->Write("gHERAH1Rho3"); gHERAH1Rho4->Write("gHERAH1Rho4"); gHERAH1Rho5->Write("gHERAH1Rho5"); hfile->Close(); cout << rootfile.c_str() <<" written." << endl; - - string parnames[16]={"B_G", "l_mass", "c_mass", "b_mass", "C", "Mu_02", "A_g", "lamda_g", "r_max", "as_qqg", "Rhs", "nLP", "R0WS", "dWS", "uCE", "EP"}; - for(int i=0; i<16; i++){ - cout<readSettingsFromFile(runcard); - settings->consolidateSettings(); + settings->readSettingsFromFile(runcard); + settings->consolidateSettings(); settings2->setA(settings->A()); settings2->setDipoleModelType(settings->dipoleModelType()); settings2->setRootfile(settings->rootfile().c_str()); - + PR(settings2->A()); PR(settings2->dipoleModelName()); PR(settings2->rootfile()); //TT - + mFlux.setS(mS); - - // if(mUseSkewednessCorrection) mbSat.useSkewednessCorrection(); + + // if(mUseSkewednessCorrection) + // mbSat.useSkewednessCorrection(); //data.SetMinQ2(0.01); //data.SetMaxQ2(651); - + //import data files //data.LoadData("/home/bharath/Project/sartre/examples/HERA_I_combined_DIS.txt", TOTAL); - data.LoadData("/home/bharath/Project/sartre/examples/HERA_II_combined_DIS.txt", TOTAL); +// data.LoadData("/Users/tbu/sartreMinuit/data/HERA_II_combined_DIS.txt", TOTAL); + data.LoadData("/star/u/ttoll/sartreMinuit/data/HERA_II_combined_DIS.txt", TOTAL); //data.LoadData("/home/bharath/Project/sartre/examples/hera_combined_sigmar_cc.txt", CHARM); MnUserParameters parameters; - //Minuit2Minimizer m2m; + //Minuit2Minimizer m2m; //Add the parameters with initial values and increment scales //////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////// MIGRAD ///////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////// //parameters.Add("B_G", params_B_G); //, 10.); //parameters.Add("light_mass", params_l_mass, 1.); //parameters.Add("charm_mass", params_c_mass, 1.); //parameters.Add("bottom_mass", params_b_mass); //parameters.Add("C", params_C, 10.); //bSat parameters //parameters.Add("mu02", params_mu02, 10.); //parameters.Add("A_g", params_A_g, 10.); //parameters.Add("lambda_g", params_lambda_g, 1.); //parameters.Add("r_max", params_r_max, 10.); //parameters.Add("as_qqg", params_as_qqg);//, 1.); //parameters.Add("Rhs", params_Rhs);//, 10.); //parameters.Add("nLP", params_nLP);//, 10.); //parameters.Add("R0WS", params_R0WS); //parameters.Add("dWS", params_dWS); //parameters.Add("uCE", params_uCE); //parameters.Add("EP", params_EP); //Vector Meson Wavefunction parameters ////parameters.Add("R_L2", RL2_jpsi);//, 10.); ////parameters.Add("N_L", NL_jpsi);//, 10.); ////parameters.Add("N_T", NT_jpsi);//, 10.); //Vector Meson enhancement parameters //parameters.Add("BV", params_BV, 10.); //parameters.Add("RV", params_RV, 10.); //parameters.Add("wV", params_wV, 10.); //parameters.Add("lambda_s", params_lambda_s); //parameters.Add("A_s", params_lambda_s); //////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////// SIMPLEX ///////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////// - parameters.Add("B_G", params_B_G); //, 0.1); - parameters.Add("light_mass", params_l_mass);//, 0.01); - parameters.Add("charm_mass", params_c_mass); //, 0.01); + parameters.Add("B_G", params_B_G);//, 0.1); + parameters.Add("light_mass", params_l_mass, 0.01); + parameters.Add("charm_mass", params_c_mass); parameters.Add("bottom_mass", params_b_mass); - parameters.Add("C", params_C); //, 0.1); + parameters.Add("C", params_C, 0.1); //bSat parameters - parameters.Add("mu02", params_mu02); //, 0.1); - parameters.Add("A_g", params_A_g); //, 0.1); - parameters.Add("lambda_g", params_lambda_g); //, 0.01); + parameters.Add("mu02", params_mu02, 0.1); + parameters.Add("A_g", params_A_g, 0.1); + parameters.Add("lambda_g", params_lambda_g, 0.01); parameters.Add("r_max", params_r_max, 0.1); parameters.Add("as_qqg", params_as_qqg); //, 0.01); parameters.Add("Rhs", params_Rhs); //, 0.5); parameters.Add("nLP", params_nLP); //, 0.1); parameters.Add("R0WS", params_R0WS); //, 1.); parameters.Add("dWS", params_dWS); //, 0.1); parameters.Add("uCE", params_uCE); //, 0.1); parameters.Add("EP", params_EP); //, 1.); - - //Vector Meson Wavefunction parameters - //parameters.Add("R_L2", RL2_jpsi); //, 0.1); - //parameters.Add("N_L", NL_jpsi); //, 0.1); - //parameters.Add("N_T", NT_jpsi); //, 0.1); - - //Vector Meson enhancement parameters - parameters.Add("BV", params_BV, 0.1); - parameters.Add("RV", params_RV, 0.1); - parameters.Add("wV", params_wV, 0.1); - - parameters.Add("lambda_s", params_lambda_s); - parameters.Add("A_s", params_lambda_s); + // + // //Vector Meson Wavefunction parameters + // //parameters.Add("R_L2", RL2_jpsi); //, 0.1); + // //parameters.Add("N_L", NL_jpsi); //, 0.1); + // //parameters.Add("N_T", NT_jpsi); //, 0.1); + // + // //Vector Meson enhancement parameters + // parameters.Add("BV", params_BV, 0.1); + // parameters.Add("RV", params_RV, 0.1); + // parameters.Add("wV", params_wV, 0.1); + // + // parameters.Add("lambda_s", params_lambda_s); + // parameters.Add("A_s", params_lambda_s); ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// //Limits for some parameters - //parameters.SetLowerLimit("A_g", 0); - //parameters.SetLowerLimit("A_s", 0); - parameters.SetLowerLimit("mu02", 0); //BS For MZ- to fix mu_02 - parameters.SetUpperLimit("mu02", 2.5); //BS For MZ- to fix mu_02 - parameters.SetLowerLimit("r_max", 0); //BS For MZ- to fix r_max - //parameters.SetUpperLimit("r_max", 15); //BS For MZ- to fix r_max +// parameters.SetUpperLimit("A_g", params_A_g); +// parameters.SetLowerLimit("A_g", params_A_g); + parameters.SetLowerLimit("mu02", 0.6);//params_mu02);//BS For MZ- to fix mu_02 + parameters.SetUpperLimit("mu02", 2.);//params_mu02); //BS For MZ- to fix mu_02 +// parameters.SetUpperLimit("r_max", params_r_max); //BS For MZ- to fix r_max +// parameters.SetLowerLimit("r_max", params_r_max); //BS For MZ- to fix r_max + parameters.SetUpperLimit("B_G", 30); parameters.SetLowerLimit("B_G", 0); - parameters.SetLowerLimit("nLP", 0); - parameters.SetLowerLimit("Rhs", 0); - // parameters.SetUpperLimit("B_G", 30); - parameters.SetLowerLimit("light_mass", 0); - parameters.SetUpperLimit("light_mass", params_c_mass); - parameters.SetLowerLimit("C", 0.); - //parameters.SetUpperLimit("C", params_C); - parameters.SetLowerLimit("A_g", 0.); - parameters.SetLowerLimit("EP", 0.); - // parameters.SetUpperLimit("EP", 1.); - parameters.SetLowerLimit("uCE", 0); - //parameters.SetUpperLimit("A_g", 0.); - //parameters.SetLowerLimit("lambda_g", params_lambda_g); - //parameters.SetUpperLimit("lambda_g", params_lambda_g); - parameters.SetLowerLimit("as_qqg", 0.); + // parameters.SetLowerLimit("nLP", 0); + //parameters.SetLowerLimit("Rhs", 0); + parameters.SetLowerLimit("light_mass", 0.);//params_l_mass); + parameters.SetUpperLimit("light_mass", 0.5);//params_l_mass); + parameters.SetLowerLimit("C", 0.5);//params_C); + parameters.SetUpperLimit("C", 6.);//params_C); + // parameters.SetLowerLimit("EP", 0.); + // // parameters.SetUpperLimit("EP", 1.); + // parameters.SetLowerLimit("uCE", 0); + // parameters.SetLowerLimit("lambda_g", params_lambda_g); + // parameters.SetUpperLimit("lambda_g", params_lambda_g); + // parameters.SetLowerLimit("as_qqg", 0.); myIntegrals.init(); DISFitter fitter(parameters); fitter.AddDataset(data); FitParameters p; p.parameter = ¶meters; vector parvec = parameters.Params(); p.values = &parvec; - + // // DGLAP Evolution can be sped up by using lookup tables // dglap.generateLookupTable(1000, 1000); dglap.useLookupTable(true); cout << "=== Initial parameters ===" << endl; cout << parameters << endl; cout << "=== Starting fit ===" << endl; // Minuit2Minimizer m2m; // cout<<"CORR is "<. # # Author: Thomas Ullrich # Last update: # $Date$ # $Author$ #=============================================================================== cmake_minimum_required (VERSION 3.1) project(SARTRE_EXAMPLES) # # Compiler options # set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long") set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_RELEASE "-O") set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) # # Set the location of the installed sartre package # set(SARTRE_DIR "$ENV{SARTRE_DIR}") if (IS_DIRECTORY "${SARTRE_DIR}" AND IS_DIRECTORY "${SARTRE_DIR}/lib" AND IS_DIRECTORY "${SARTRE_DIR}/include" AND IS_DIRECTORY "${SARTRE_DIR}/gemini" AND IS_DIRECTORY "${SARTRE_DIR}/tables") message (STATUS "Using Sartre installed in ${SARTRE_DIR}") else() message(FATAL_ERROR "Environment variable SARTRE_DIR does not point to a valid Sartre installation.") endif() message(STATUS "System = ${CMAKE_SYSTEM}") # # User can select if Sartre is to be compiled # in multithreaded mode or not. # Default is off. # cmake -DMULTITHREADED=ON . # option(MULTITHREADED "Switch on/off multithreaded mode" OFF) if (MULTITHREADED) add_definitions(-DSARTRE_IN_MULTITHREADED_MODE) message (STATUS "MULTITHREADED mode is ON") endif (MULTITHREADED) # # See if optional libraries (Glut & OpenGL) are present. # If not skip the examples that require them. # # GLUT #find_package(GLUT) #if (GLUT_FOUND) # message("blah") # include_directories(${GLUT_INCLUDE_DIR}) # set(LIBS ${LIBS} ${GLUT_LIBRARIES}) #endif() # OpenGL find_package(OpenGL) if (OPENGL_FOUND) include_directories(${OPENGL_INCLUDE_DIR}) set(LIBS ${LIBS} ${OPENGL_LIBRARIES}) endif() # # Built all programs in example directory # add_executable(createBSatBDependenceTable createBSatBDependenceTable.cpp) add_executable(heraCompare heraCompare.cpp) add_executable(sartreMain sartreMain.cpp) add_executable(sartreTest sartreTest.cpp) add_executable(tableGeneratorMain tableGeneratorMain.cpp) add_executable(tableGeneratorUPCMain tableGeneratorUPCMain.cpp) add_executable(lambdaGeneratorMain lambdaGeneratorMain.cpp) add_executable(sartreMinuitMain sartreMinuitMain.cpp) #if (OPENGL_FOUND AND GLUT_FOUND) # add_executable(nucleusAnimation nucleusAnimation.cpp) #endif() # # Make sure the files to search for external libraries are available # set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/modules") # # Compiler options # set(CMAKE_C_FLAGS_DEBUG " -g -W") set(CMAKE_CXX_FLAGS_DEBUG " -g -W -Wall -Wextra -pedantic -Wno-long-long -std=c++11") set(CMAKE_C_FLAGS_RELEASE " -O -W") set(CMAKE_CXX_FLAGS_RELEASE " -O -W -Wall -Wextra -pedantic -Wno-long-long -std=c++11") set(CMAKE_C_FLAGS " -O -W") set(CMAKE_CXX_FLAGS " -O -std=c++11 -W -Wall -Wextra -pedantic -Wno-long-long") # # Include files locations # include_directories("${SARTRE_DIR}/include") # # Find external required packages # (see also FindGSL.cmake and FindROOT.cmke in cmake/modules) # # GSL find_package(GSL REQUIRED) include_directories(${GSL_INCLUDE_DIR}) set(LIBS ${LIBS} ${GSL_LIBRARIES}) # ROOT find_package(ROOT REQUIRED) include_directories(${ROOT_INCLUDE_DIR}) -set(LIBS ${LIBS} ${ROOT_LIBRARIES} "-lEG -lMathMore -lUnuran" ) +set(LIBS ${LIBS} ${ROOT_LIBRARIES} "-lEG -lMathMore -lMinuit2 -lUnuran" ) #BOOST if (MULTITHREADED) set(Boost_USE_MULTITHREADED ON) find_package(Boost 1.39 COMPONENTS thread REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIR}) - set(LIBS ${LIBS} ${Boost_LIBRARIES}) + set(LIBS ${LIBS} ${Boost_LIBRARIES} -lboost_system) endif(Boost_FOUND) endif (MULTITHREADED) # # Libraries to link with # #link_directories(${SARTRE_DIR}/lib) file (GLOB ALLSARTRE ${SARTRE_DIR}/lib/*) set(LIBS ${ALLSARTRE} ${ALLSARTRE} ${LIBS}) target_link_libraries(createBSatBDependenceTable ${LIBS}) target_link_libraries(heraCompare ${LIBS}) target_link_libraries(sartreMain ${LIBS}) target_link_libraries(sartreTest ${LIBS}) target_link_libraries(tableGeneratorMain ${LIBS}) target_link_libraries(tableGeneratorUPCMain ${LIBS}) target_link_libraries(lambdaGeneratorMain ${LIBS}) target_link_libraries(sartreMinuitMain ${LIBS}) #if (OPENGL_FOUND AND GLUT_FOUND) # target_link_libraries(nucleusAnimation ${LIBS}) #endif()