Page MenuHomeHEPForge
Files F8308822

No OneTemporary
Actions

View Options

Index: appl-bridge/trunk/src/standpdfCorr.cxx
===================================================================
--- appl-bridge/trunk/src/standpdfCorr.cxx (revision 1322)
+++ appl-bridge/trunk/src/standpdfCorr.cxx (revision 1323)
@@ -1,326 +1,327 @@
#include <iostream>
#include <string>
#include <vector>
#include <stdio.h>
#include <cstdlib>
#include "appl_grid/appl_grid.h"
#include "appl_grid/Directory.h"
#include "appl_grid/appl_timer.h"
#include <TCanvas.h>
#include <TString.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TPad.h>
#include <TString.h>
#include <TMath.h>
#include "LHAPDF/LHAPDF.h"
#include "hoppet_v1.h"
#define DBG false
static const int nFlavours = 5;
static const Int_t nCorrPDFs = 5;
static const TString corrPDFName[nCorrPDFs] = { "gluon" , "d_val" , "u_val", "strange", "sea"};
// static const int nPDFs = 1;
// static const TString PDFName[nPDFs] =
// {
// "NNPDF11_100"
// };
// static const int nPDFMembers[nPDFs]={100};
// static const int nPDFMembersMax = 1000;
static const int nPDFs = 1;
static const TString PDFName[nPDFs] =
{
// "NNPDF21_100",
"NNPDF23_nlo_as_0119",
// "NNPDF23_nlo_noLHC_as_0118",
};
//static const int nPDFMembers[nPDFs] = { 101, 101, 41, 21, 27, 29 };
static const int nPDFMembers[nPDFs] = { 100 };
static const int nPDFMembersMax = 100;
// wrapper to get the x*pdf from gavins evolution code.
// in fact, don't actually need this wrapper any longer - we could
// just pass the routine directly
extern "C" void evolvepdf_(const double& x, const double& Q, double* xf);
extern "C" double alphaspdf_(const double& Q);
void GetPdf(const double& x, const double& Q, double* xf) {
double _x_(x);
if (_x_ < 1e-7)_x_ = 1e-7;
evolvepdf_( _x_, Q, xf);
xf[0] = 0; xf[12] = 0;
return;
}
// void GetPdfSplit(const double& x, const double& Q, double* f)
// {
// hoppetevalsplit_( x, Q, nLoops, nFlavours, f);
// return;
// }
TH1D* divide( const TH1D* h1, const TH1D* h2 ) {
if ( h1==NULL || h2==NULL ) return NULL;
TH1D* h = (TH1D*)h1->Clone();
if ( DBG ) std::cout << "histograms " << h1->GetTitle() << " " << h2->GetTitle() << std::endl;
for ( int i=1 ; i<=h1->GetNbinsX() ; i++ ) {
double b = h2->GetBinContent(i);
//double be = h2->GetBinError(i);
double t = h1->GetBinContent(i);
//double te = h1->GetBinError(i);
double r = ( b!=0 ? t/b : 0 );
// double re = ( b!=0 ? sqrt((r+1)*r/b) : 0 );
double re = 0;
h->SetBinContent( i, r );
h->SetBinError( i, re ) ;
// if ( DBG ) std::cout << "\tx=" << h->GetBinCenter(i) << "\tratio=" << r << std::endl;
}
double hmin = h->GetBinContent(1);
double hmax = h->GetBinContent(1);
for ( int i=2 ; i<=h->GetNbinsX() ; i++ ) {
double d = h->GetBinContent(i);
if ( hmin>d ) hmin=d;
if ( hmax<d ) hmax=d;
}
if ( DBG ) std::cout << "\tmin ratio = " << hmin << "\tmax ratio = " << hmax << std::endl;
if ( h->GetMaximum()<1.01 ) h->SetMaximum(1.01);
if ( h->GetMinimum()>0.99 ) h->SetMinimum(0.99);
return h;
}
Double_t getSigma(std::vector<double>& v , Double_t& val)
{
Double_t sigma(0);
for(UInt_t i=0; i< v.size(); i++)
sigma += TMath::Power( val - v[i],2)/ v.size();
sigma = TMath::Sqrt( sigma );
return sigma;
}
Double_t getAverage(std::vector<double> &v )
{
Double_t ave(0);
for(UInt_t i=0; i< v.size(); i++) ave+= v[i]/v.size();
return ave;
}
// extern double Escale;
int main(int argc, char** argv) {
if ( argc<2 ) return -1;
TString gridName(argv[1]);
appl::grid g(gridName.Data());
g.trim();
const int nLoops = g.nloops();
+// const int nLoops = 1;
// get grid cms energy
double gridEnergy = g.getCMSScale();
//double Escale = 1.;
std::cout << "grid CMS energy " << gridEnergy << std::endl;
// get all the reference histograms
TFile* f;
f = new TFile(argv[1]);
TString outFileName = gridName(0,gridName.Index(".root"));
outFileName += "_standPDFCorr.root";
TFile* fout = new TFile(outFileName, "recreate");
TH1D* reference = (TH1D*)f->Get("grid/reference");
const Int_t nxBins = 100;
Double_t xBins[nxBins+1];
for (int iBin=1; iBin<=nxBins;iBin++) xBins[iBin]=double(iBin)/double(nxBins);
const Int_t nq2Bins = reference->GetNbinsX();
Double_t q2Bins[nq2Bins+1];
for (int iBin=0; iBin<=nq2Bins;iBin++) q2Bins[iBin]=reference->GetBinLowEdge(iBin+1);
TString histName("");
TH1D *xsAverage, *xsReplica[nPDFMembersMax];
TH2D* rhoHist[nCorrPDFs], *actualPDFsXQ[13][nPDFMembersMax];
Double_t actualPDFs[13];
std::vector<double> sigma, pdfValues, sigmaTimesPDF;
sigma. reserve(nPDFMembersMax);
pdfValues. reserve(nPDFMembersMax);
sigmaTimesPDF.reserve(nPDFMembersMax);
// now calculate all the cross sections
for (int iPDF = 0; iPDF < nPDFs; iPDF++)
{
TString PDFSetName = "PDFsets/";
PDFSetName += PDFName[iPDF];
PDFSetName += ".LHgrid";
LHAPDF::initPDFSet(PDFSetName.Data(), 0);
Directory pdfDir(PDFName[iPDF].Data());
pdfDir.push();
xsAverage = g.convolute(GetPdf, alphaspdf_ , nLoops);
std::cout << "\t\t average done .... " << std::endl;
histName.Form("xs_Average_%s", PDFName[iPDF].Data() );
xsAverage->SetName( histName );
xsAverage->SetTitle( histName );
xsAverage->Write("", TObject::kOverwrite);
Directory pdfmemberDir("replicas");
pdfmemberDir.push();
for ( int iReplica = 0; iReplica < nPDFMembers[iPDF]; iReplica++ )
{
std::cout << " working with replica "<< iReplica << std::endl;
LHAPDF::usePDFMember( iReplica + 1 );
xsReplica[iReplica] = g.convolute(GetPdf, alphaspdf_ , nLoops );
histName.Form("xs_Average_%s_replica_%d", PDFName[iPDF].Data(), iReplica );
xsReplica[iReplica] -> SetName( histName );
xsReplica[iReplica] -> SetTitle( histName );
xsReplica[iReplica] -> Write("", TObject::kOverwrite);
for (int j=0; j<13; j++)
{
histName.Form("actualPDF_flavor_%d_replica_%d", j, iReplica);
actualPDFsXQ[j][iReplica] = new TH2D(histName, histName, nxBins, xBins, nq2Bins, q2Bins);
actualPDFsXQ[j][iReplica]->SetDirectory(0);
}
for (int ixbin=0; ixbin< nxBins; ixbin++)
for (int iqbin=0; iqbin< nq2Bins; iqbin++)
{
//for (int ii=0; ii<13;ii++) actualPDFs[ii] = 0;
double x = (xBins[ixbin+1] + xBins[ixbin])/2;
double q = (q2Bins[iqbin+1] + q2Bins[iqbin])/2;
GetPdf( x, q, actualPDFs );
for (int j=0; j<13; j++) actualPDFsXQ[j][iReplica]->SetBinContent( ixbin+1, iqbin+1, actualPDFs[j] );
}
}
pdfmemberDir.pop();
//for (int j=0; j<13; j++) actualPDFsXQ[j]->Write("", TObject::kOverwrite);
std::cout << "\t\t pdfs done .... " << std::endl;
std::cout << "\t\t starting correlations .... " << std::endl;
Directory correlationDir("correlations");
correlationDir.push();
for (int iCorr = 0; iCorr < nCorrPDFs; iCorr++)
{
histName.Form("correlation_%s", corrPDFName[iCorr].Data() );
rhoHist[iCorr] = new TH2D(histName, histName, nxBins, xBins, nq2Bins, q2Bins);
for (int ixbin=0; ixbin< nxBins; ixbin++)
for (int iqbin=0; iqbin< nq2Bins; iqbin++)
{
sigma.clear();
pdfValues.clear();
sigmaTimesPDF.clear();
for ( int iReplica = 0; iReplica < nPDFMembers[iPDF]; iReplica++ )
{
Double_t valS = xsReplica[iReplica]->GetBinContent(iqbin+1);
Double_t valP = 1e6;
if (iCorr == 0)
valP = actualPDFsXQ[6][iReplica]->GetBinContent( ixbin+1, iqbin+1);
else if (iCorr == 1)
valP = (
actualPDFsXQ[7][iReplica]->GetBinContent( ixbin+1, iqbin+1)-
actualPDFsXQ[5][iReplica]->GetBinContent( ixbin+1, iqbin+1)
);
else if (iCorr == 2)
valP = (
actualPDFsXQ[8][iReplica]->GetBinContent( ixbin+1, iqbin+1)-
actualPDFsXQ[4][iReplica]->GetBinContent( ixbin+1, iqbin+1)
);
else if (iCorr == 3)
valP = actualPDFsXQ[9][iReplica]->GetBinContent( ixbin+1, iqbin+1);
else if (iCorr == 4 )
valP = (
actualPDFsXQ[3][iReplica]->GetBinContent( ixbin+1, iqbin+1)+
actualPDFsXQ[4][iReplica]->GetBinContent( ixbin+1, iqbin+1)+
actualPDFsXQ[5][iReplica]->GetBinContent( ixbin+1, iqbin+1)
);
else
{
std::cout <<"Unknown iCorr !!! "<< std::endl;
exit(-1);
}
sigma.push_back( valS );
pdfValues.push_back( valP );
sigmaTimesPDF.push_back( valS * valP );
}
Double_t aveS = getAverage(sigma);
Double_t aveP = getAverage(pdfValues);
Double_t rhoVal = getAverage(sigmaTimesPDF) - aveS*aveP;
rhoVal /= getSigma(sigma, aveS)*getSigma(pdfValues, aveP);
rhoVal *= double(nPDFMembers[iPDF])/double(nPDFMembers[iPDF]-1);
rhoHist[iCorr] -> SetBinContent( ixbin+1, iqbin+1, rhoVal );
}
rhoHist[iCorr]->SetDrawOption("colz");
rhoHist[iCorr]->Write("", TObject::kOverwrite);
}
correlationDir.pop();
std::cout << "\t\t correlations done.... " << std::endl;
pdfDir.pop();
}
fout->Write();
fout->Close();
return 0;
}
Index: appl-bridge/trunk/src/standpdf.cxx
===================================================================
--- appl-bridge/trunk/src/standpdf.cxx (revision 1322)
+++ appl-bridge/trunk/src/standpdf.cxx (revision 1323)
@@ -1,431 +1,432 @@
#include <iostream>
#include <string>
#include <vector>
#include <stdio.h>
#include <cstdlib>
#include "appl_grid/appl_grid.h"
#include "appl_grid/Directory.h"
#include "appl_grid/appl_timer.h"
#include <TCanvas.h>
#include <TString.h>
#include <TH1D.h>
#include <TFile.h>
#include <TPad.h>
#include "LHAPDF/LHAPDF.h"
#include "hoppet_v1.h"
#define DBG false
static const int nFlavours = 5;
const int Nscales=7;
const int nScales=Nscales;
static const double mur[Nscales] = { 1.0, 0.5, 2.0, 1.0, 1.0, 0.5, 2.0 };
static const double muf[Nscales] = { 1.0, 1.0, 1.0, 0.5, 2.0, 0.5, 2.0 };
// static const int nPDFs = 1;
// static const TString PDFName[nPDFs] =
// {
// "NNPDF11_100"
// };
// static const int nPDFMembers[nPDFs]={100};
// static const int nPDFMembersMax = 1000;
static const int nPDFs = 9;
static const TString PDFName[nPDFs] =
{
"cteq66",
"CT10nlo",
"NNPDF21_100",
"NNPDF23_nlo_as_0119",
"NNPDF23_nlo_noLHC_as_0118",
"MSTW2008nlo90cl",
"HERAPDF15NLO_EIG",
"GJR08VFnloE",
"abm11_5n_nlo"
};
//static const int nPDFMembers[nPDFs] = { 101, 101, 41, 21, 27, 29 };
static const int nPDFMembers[nPDFs] = {45, 53, 101, 101, 101, 41, 21, 27, 29 };
static const int nPDFMembersMax = 101;
extern "C" void evolvepdf_(const double& x, const double& Q, double* xf);
extern "C" double alphaspdf_(const double& Q);
void GetPdf(const double& x, const double& Q, double* xf) {
double _x_(x);
if (_x_ < 1e-7)_x_ = 1e-7;
evolvepdf_( _x_, Q, xf);
xf[0] = 0; xf[12] = 0;
return;
}
TH1D* divide( const TH1D* h1, const TH1D* h2 ) {
if ( h1==NULL || h2==NULL ) return NULL;
TH1D* h = (TH1D*)h1->Clone();
if ( DBG ) std::cout << "histograms " << h1->GetTitle() << " " << h2->GetTitle() << std::endl;
for ( int i=1 ; i<=h1->GetNbinsX() ; i++ ) {
double b = h2->GetBinContent(i);
//double be = h2->GetBinError(i);
double t = h1->GetBinContent(i);
//double te = h1->GetBinError(i);
double r = ( b!=0 ? t/b : 0 );
// double re = ( b!=0 ? sqrt((r+1)*r/b) : 0 );
double re = 0;
h->SetBinContent( i, r );
h->SetBinError( i, re ) ;
// if ( DBG ) std::cout << "\tx=" << h->GetBinCenter(i) << "\tratio=" << r << std::endl;
}
double hmin = h->GetBinContent(1);
double hmax = h->GetBinContent(1);
for ( int i=2 ; i<=h->GetNbinsX() ; i++ ) {
double d = h->GetBinContent(i);
if ( hmin>d ) hmin=d;
if ( hmax<d ) hmax=d;
}
if ( DBG ) std::cout << "\tmin ratio = " << hmin << "\tmax ratio = " << hmax << std::endl;
if ( h->GetMaximum()<1.01 ) h->SetMaximum(1.01);
if ( h->GetMinimum()>0.99 ) h->SetMinimum(0.99);
return h;
}
// extern double Escale;
int main(int argc, char** argv) {
if ( argc<2 ) return -1;
TString gridName(argv[1]);
appl::grid g(gridName.Data());
g.trim();
const int nLoops = g.nloops();
+// const int nLoops = 1;
// get grid cms energy
double gridEnergy = g.getCMSScale();
//double Escale = 1.;
std::cout << "grid CMS energy " << gridEnergy << std::endl;
// get all the reference histograms
//TFile* f;
//f = new TFile(argv[1]);
TString outFileName = gridName(0,gridName.Index(".root"));
outFileName += "_standTheoryUncert.root";
TFile* fout = new TFile(outFileName, "recreate");
// now calculate all the cross sections
for (int iPDF = 0; iPDF < nPDFs; iPDF++)
{
TString PDFSetName = "PDFsets/";
PDFSetName += PDFName[iPDF];
PDFSetName += ".LHgrid";
LHAPDF::initPDFSet(PDFSetName.Data(), 0);
Directory pdfDir(PDFName[iPDF].Data());
pdfDir.push();
Directory scalesDir("scales");
scalesDir.push();
TH1D* htemp = NULL;
TH1D* xsec_scale[Nscales];
TH1D* ratio_scale[Nscales];
for ( int i = 0; i < Nscales ; i++ )
{
char hname[64];
sprintf( hname, "xsec_scale_%d", i);
htemp = g.convolute(GetPdf, alphaspdf_ , nLoops, mur[i], muf[i]);
xsec_scale[i] = htemp;
//xsec_scale[i] = (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_scale[i]->SetName(hname);
xsec_scale[i]->Write();
sprintf( hname, "ratio_scale_%d", i);
ratio_scale[i] = divide( xsec_scale[i], xsec_scale[0] );
if ( ratio_scale[i] )
{
ratio_scale[i]->SetName(hname);
ratio_scale[i]->Write("", TObject::kOverwrite);
}
}
scalesDir.pop();
std::cout << "\t\t scales done .... " << std::endl;
// now take all the ratios etc
Directory pdfmemberDir("pdfs");
pdfmemberDir.push();
TH1D* xsec_pdf[nPDFMembersMax];
TH1D* ratio_pdf[nPDFMembersMax];
for ( int i = 0; i < nPDFMembers[iPDF]; i++ )
{
LHAPDF::usePDFMember(i);
char hname[64];
sprintf( hname, "xsec_pdf_%d", i);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_pdf[i] = htemp;
//xsec_pdf[i] = (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_pdf[i]->SetName(hname);
xsec_pdf[i]->Write("", TObject::kOverwrite);
sprintf( hname, "ratio_pdf_%d", i);
ratio_pdf[i] = divide( xsec_pdf[i], xsec_pdf[0]);
if ( ratio_pdf[i] )
{
ratio_pdf[i]->SetName(hname);
ratio_pdf[i]->Write("", TObject::kOverwrite);
}
}
pdfmemberDir.pop();
std::cout << "\t\t pdfs done .... " << std::endl;
//
// alpha2 uncertainty
//
TString xsalpha_Name[2] = {"xsec_alphaMinus", "xsec_alphaPlus"};
TString ratioalpha_Name[2] = {"ratio_alphaMinus", "ratio_alphaPlus"};
TH1D* xsec_alphas[2];
TH1D* ratio_alphas[2];
Directory alphasDir("alphas");
alphasDir.push();
if (iPDF == 0) // cteq66
{
LHAPDF::initPDFSet("PDFsets/cteq66alphas.LHgrid",4);
LHAPDF::usePDFMember(4); // plus 0.120
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1] = htemp;
//xsec_alphas [1] = (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::usePDFMember(0); //minus 0.116
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 1) // ct10
{
LHAPDF::initPDFSet("PDFsets/CT10nlo_as_0120.LHgrid",0);
LHAPDF::usePDFMember(0); // plus 0.120
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1] = htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::initPDFSet("PDFsets/CT10nlo_as_0116.LHgrid",0);
LHAPDF::usePDFMember(0); //minus 0.116
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 2) // nnpdf21
{
LHAPDF::initPDFSet("PDFsets/NNPDF21_as_0121_100.LHgrid",0); //plus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1]= htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::initPDFSet("PDFsets/NNPDF21_as_0117_100.LHgrid",0); //minus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 3) // nnpdf23
{
LHAPDF::initPDFSet("PDFsets/NNPDF23_nlo_as_0121.LHgrid",0); //plus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1]= htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::initPDFSet("PDFsets/NNPDF23_nlo_as_0117.LHgrid",0); //minus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 4) // nnpdf23_noLHC
{
LHAPDF::initPDFSet("PDFsets/NNPDF23_nlo_noLHC_as_0120.LHgrid",0); //plus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1]= htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::initPDFSet("PDFsets/NNPDF23_nlo_noLHC_as_0116.LHgrid",0); //minus
LHAPDF::usePDFMember(0);
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 5) // mstw
{
LHAPDF::initPDFSet("PDFsets/MSTW2008nlo90cl_asmz+90cl.LHgrid",0);
LHAPDF::usePDFMember(0); //plus
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1]= htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::initPDFSet("PDFsets/MSTW2008nlo90cl_asmz-90cl.LHgrid",0);
LHAPDF::usePDFMember(0); //minus
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else if (iPDF == 6) // hera
{
LHAPDF::initPDFSet("PDFsets/HERAPDF15NLO_ALPHAS.LHgrid",11);
LHAPDF::usePDFMember(11); //plus
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [1]= htemp;
//xsec_alphas [1]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [1]->SetName(xsalpha_Name[1]);
ratio_alphas[1] = divide( xsec_alphas[1], xsec_pdf[0]);
ratio_alphas[1]->SetName(ratioalpha_Name[1]);
//----------------
LHAPDF::usePDFMember(9); //minus
htemp = g.convolute(GetPdf, alphaspdf_);
xsec_alphas [0]= htemp;
//xsec_alphas [0]= (TH1D*)htemp->Rebin(_nRebins, TString::Format("%s_rebin",htemp->GetName()), _binsRebin);
//delete htemp;
xsec_alphas [0]->SetName(xsalpha_Name[0]);
ratio_alphas[0] = divide( xsec_alphas[0], xsec_pdf[0]);
ratio_alphas[0]->SetName(ratioalpha_Name[0]);
}
else
{
xsec_alphas [0] = (TH1D*)xsec_pdf[0]->Clone(xsalpha_Name[0]);
xsec_alphas [1] = (TH1D*)xsec_pdf[0]->Clone(xsalpha_Name[1]);
ratio_alphas [0] = (TH1D*)xsec_pdf[0]->Clone(ratioalpha_Name[0]);
ratio_alphas [1] = (TH1D*)xsec_pdf[0]->Clone(ratioalpha_Name[1]);
for (int i = 0; i < 2; i++)
{
xsec_alphas [i]->Reset();
ratio_alphas [i]->Reset();
}
}
for (int i = 0; i < 2; i++)
{
if (xsec_alphas[i])
{
xsec_alphas[i] -> Write("", TObject::kOverwrite);
delete xsec_alphas[i];
}
if (ratio_alphas[i])
{
ratio_alphas[i] -> Write("", TObject::kOverwrite);
delete ratio_alphas[i];
}
}
for ( int i = 0; i < nPDFMembers[iPDF]; i++ )
{
if (ratio_pdf[i]) delete ratio_pdf[i];
if (xsec_pdf[i]) delete xsec_pdf[i];
}
for ( int i = 0; i < Nscales ; i++ )
{
if (ratio_scale[i]) delete ratio_scale[i];
if(xsec_scale[i]) delete xsec_scale[i];
}
alphasDir.pop();
pdfDir.pop();
std::cout << "\t\t alphas done .... " << std::endl;
}
fout->Write();
fout->Close();
return 0;
}
Index: appl-bridge/trunk/src/standsubproc.cxx
===================================================================
--- appl-bridge/trunk/src/standsubproc.cxx (revision 1322)
+++ appl-bridge/trunk/src/standsubproc.cxx (revision 1323)
@@ -1,282 +1,318 @@
#include <iostream>
#include <string>
#include <vector>
#include <stdio.h>
#include <cstdlib>
#include "appl_grid/appl_grid.h"
#include "appl_grid/Directory.h"
#include "appl_grid/appl_timer.h"
#include <TCanvas.h>
#include <TH1D.h>
#include <TFile.h>
#include <TPad.h>
#include "LHAPDF/LHAPDF.h"
#include "hoppet_v1.h"
#define DBG false
static const int nFlavours = 5;
extern "C" void evolvepdf_(const double& x, const double& Q, double* xf);
extern "C" double alphaspdf_(const double& Q);
// wrapper to get the x*pdf from gavins evolution code.
// in fact, don't actually need this wrapper any longer - we could
// just pass the routine directly
void GetPdf(const double& x, const double& Q, double* xf) {
double _x_(x);
if (_x_ < 1e-7)_x_ = 1e-7;
evolvepdf_( _x_, Q, xf);
xf[0] = 0; xf[12] = 0;
return;
}
static int iFlav(-1000), jFlav(-1000);
void GetPdfP1(const double& x, const double& Q, double* xf) {
double _x_(x);
if (_x_ < 1e-7)_x_ = 1e-7;
evolvepdf_( _x_, Q, xf);
xf[0] = 0; xf[12] = 0;
if ( std::abs(iFlav) > nFlavours ) return;
for (int f = -6; f <= 6; f++) if (f != iFlav) xf[f+6] = 0;
return;
}
void GetPdfP2(const double& x, const double& Q, double* xf) {
double _x_(x);
if (_x_ < 1e-7)_x_ = 1e-7;
evolvepdf_( _x_, Q, xf);
xf[0] = 0; xf[12] = 0;
if ( std::abs(jFlav) > nFlavours ) return;
for (int f = -6; f <= 6; f++) if (f != jFlav) xf[f+6] = 0;
return;
}
// void GetPdfSplit(const double& x, const double& Q, double* f)
// {
// hoppetevalsplit_( x, Q, nLoops, nFlavours, f);
// return;
// }
void increment( TH1D* h1, const TH1D* h2 ) {
if ( h1==NULL || h2==NULL ) return;
if ( DBG ) std::cout << "increment histograms " << h1->GetTitle() << " " << h2->GetTitle() << std::endl;
for ( int i=1 ; i<=h1->GetNbinsX() ; i++ ) {
double b1 = h1->GetBinContent(i);
double b2 = h2->GetBinContent(i);
h1->SetBinContent(i, b1+b2);
}
}
TH1D* divide( const TH1D* h1, const TH1D* h2 ) {
if ( h1==NULL || h2==NULL ) return NULL;
TH1D* h = (TH1D*)h1->Clone();
if ( DBG ) std::cout << "histograms " << h1->GetTitle() << " " << h2->GetTitle() << std::endl;
for ( int i=1 ; i<=h1->GetNbinsX() ; i++ ) {
double b = h2->GetBinContent(i);
//double be = h2->GetBinError(i);
double t = h1->GetBinContent(i);
//double te = h1->GetBinError(i);
double r = ( b!=0 ? t/b : 0 );
// double re = ( b!=0 ? sqrt((r+1)*r/b) : 0 );
double re = 0;
h->SetBinContent( i, r );
h->SetBinError( i, re ) ;
// if ( DBG ) std::cout << "\tx=" << h->GetBinCenter(i) << "\tratio=" << r << std::endl;
}
double hmin = h->GetBinContent(1);
double hmax = h->GetBinContent(1);
for ( int i=2 ; i<=h->GetNbinsX() ; i++ ) {
double d = h->GetBinContent(i);
if ( hmin>d ) hmin=d;
if ( hmax<d ) hmax=d;
}
if ( DBG ) std::cout << "\tmin ratio = " << hmin << "\tmax ratio = " << hmax << std::endl;
if ( h->GetMaximum()<1.01 ) h->SetMaximum(1.01);
if ( h->GetMinimum()>0.99 ) h->SetMinimum(0.99);
return h;
}
-const int nPDFs = 1;
+const int nPDFs = 10;
const TString pdfNames [nPDFs] =
{
"CT10nlo"
- // , "cteq66", "NNPDF21_100",
- // "MSTW2008lo68cl",
- // "MSTW2008nlo68cl",
- // "MSTW2008nnlo68cl"
- //, "MSTW2008nlo_NF3", "MSTW2008nlo_NF4","MRST2007lomod", "MRSTMCal"
+ , "cteq66", "NNPDF21_100",
+ "MSTW2008lo68cl",
+ "MSTW2008nlo68cl",
+ "MSTW2008nnlo68cl"
+ , "MSTW2008nlo_NF3", "MSTW2008nlo_NF4","MRST2007lomod", "MRSTMCal"
};
+static const int nPDFMembers[nPDFs] = { 53, 45, 101, 41, 41, 41, 41, 41, 0, 0 };
+static const int nPDFMembersMax = 101;
// const TString pdfNamesFull [nPDFs] =
// {
// "PDFsets/CT10nlo.LHgrid", "PDFsets/cteq66.LHgrid", "PDFsets/NNPDF21_100.LHgrid",
// "PDFsets/MSTW2008nlo68cl.LHgrid", "PDFsets/MSTW2008nlo90cl_nf3.LHgrid", "PDFsets/MSTW2008nlo90cl_nf4.LHgrid",
// "PDFsets/MRST2007lomod.LHgrid", "PDFsets/MRSTMCal.LHgrid"
// };
// extern double Escale;
int main(int argc, char** argv) {
if ( argc<2 ) return -1;
TString gridName(argv[1]);
appl::grid g(gridName.Data());
g.trim();
const int nLoops = g.nloops();
+// const int nLoops = 1;
// get grid cms energy
double gridEnergy = g.getCMSScale();
std::cout << "grid CMS energy " << gridEnergy << std::endl;
// get all the reference histograms
TString outFileName = gridName(0,gridName.Index(".root"));
outFileName += "_standSubproc.root";
TFile* fout = new TFile(outFileName, "recreate");
const int nSubProcesses = g.subProcesses(nLoops);
iFlav=-1000; jFlav=-1000;
for (int iPDF = 0; iPDF < nPDFs; iPDF++)
{
TString pdfNameFull("");
pdfNameFull.Form("PDFsets/%s.LHgrid",pdfNames[iPDF].Data());
LHAPDF::initPDFSet( pdfNameFull.Data(), 0 );
Int_t _nLoops_(nLoops);
//if ( (iPDF == 6) || (iPDF == 7) ) _nLoops_ = 0;
Directory pdfDir(pdfNames[iPDF].Data());
pdfDir.push();
TH1D* xsec = g.convolute(GetPdf, alphaspdf_, _nLoops_);
//std::cout << "simple conv " << std::endl;
xsec->SetName("xsec");
xsec->SetTitle("xsec");
xsec->Write("",TObject::kOverwrite);
TH1D* test = (TH1D*)xsec->Clone("testSubProcesses");
test->Reset();
TH1D* testp = (TH1D*)xsec->Clone("testSubProcessesPartons");
testp->Reset();
Directory subDir("subProcesses");
subDir.push();
TH1D* xsecsub[nSubProcesses];
TH1D* xsecsubRatio[nSubProcesses];
for (int isp = 0; isp < nSubProcesses; isp++)
{
xsecsub[isp] = g.convolute_subproc(isp, GetPdf, alphaspdf_, _nLoops_);
TString xsecsubName("");
xsecsubName.Form("xsec_subproc_%03d",isp);
xsecsub[isp]->SetName(xsecsubName);
xsecsub[isp]->SetTitle(xsecsubName);
test->Add(xsecsub[isp]);
xsecsubRatio[isp] = divide(xsecsub[isp], xsec);
xsecsubName.Form("xsec_subproc_%03d_ratio",isp);
xsecsubRatio[isp]->SetName(xsecsubName);
xsecsubRatio[isp]->SetTitle(xsecsubName);
xsecsub[isp]->Write("",TObject::kOverwrite);
xsecsubRatio[isp]->Write("",TObject::kOverwrite);
//std::cout << "subproc conv " <<isp<< std::endl;
+ delete xsecsub[isp];
+ delete xsecsubRatio[isp];
+
}
+ //delete [] xsecsub;
+ //delete [] xsecsubRatio;
subDir.pop();
test->Divide(xsec);
test->Write("",TObject::kOverwrite);
-
+ delete test;
Directory subpDir("subProcessesPartons");
subpDir.push();
TH1D* xsecsubP[2*nFlavours + 1][2*nFlavours + 1];
TH1D* xsecsubPRatio[2*nFlavours + 1][2*nFlavours + 1];
-
+
for (int i = -nFlavours; i <= nFlavours; i++)
for (int j = -nFlavours; j <= nFlavours; j++)
{
- iFlav=i; jFlav=j;
+ LHAPDF::usePDFMember(0);
+
+ iFlav=i; jFlav=j;
xsecsubP[i][j] = g.convolute(GetPdfP1, GetPdfP2, alphaspdf_, _nLoops_);
TString xsecsubName("");
xsecsubName.Form("xsec_subprocPartons_%02d_%02d",i,j);
xsecsubP[i][j]->SetName(xsecsubName);
xsecsubP[i][j]->SetTitle(xsecsubName);
testp->Add(xsecsubP[i][j]);
xsecsubPRatio[i][j] = divide(xsecsubP[i][j], xsec);
xsecsubName.Form("xsec_subprocPartons_%02d_%02d_ratio",i,j);
xsecsubPRatio[i][j]->SetName(xsecsubName);
xsecsubPRatio[i][j]->SetTitle(xsecsubName);
xsecsubP[i][j]->Write("",TObject::kOverwrite);
xsecsubPRatio[i][j]->Write("",TObject::kOverwrite);
- //std::cout << "subproc conv2 " <<i<<" "<<j<< std::endl;
+ delete xsecsubP[i][j];
+ delete xsecsubPRatio[i][j];
+
+ xsecsubName.Form("PDFUncertainties_%02d_%02d",i,j);
+ Directory subpuncertDir(xsecsubName.Data());
+ subpuncertDir.push();
+ for ( int iMem = 1; iMem < nPDFMembers[iPDF]; iMem++ )
+ {
+ LHAPDF::usePDFMember( iMem );
+ xsecsubP[i][j] = g.convolute(GetPdfP1, GetPdfP2, alphaspdf_, _nLoops_);
+ xsecsubName.Form("xsec_subprocPartons_%02d_%02d_var_%03d",i,j,iMem);
+ xsecsubP[i][j]->SetName(xsecsubName);
+ xsecsubP[i][j]->SetTitle(xsecsubName);
+
+ xsecsubPRatio[i][j] = divide(xsecsubP[i][j], xsec);
+ xsecsubName.Form("xsec_subprocPartons_%02d_%02d_var_%03d_ratio",i,j,iMem);
+ xsecsubPRatio[i][j]->SetName(xsecsubName);
+ xsecsubPRatio[i][j]->SetTitle(xsecsubName);
+
+ xsecsubP[i][j]->Write("",TObject::kOverwrite);
+ xsecsubPRatio[i][j]->Write("",TObject::kOverwrite);
+ delete xsecsubP[i][j];
+ delete xsecsubPRatio[i][j];
+ }
+
+ subpuncertDir.pop();
+
}
subpDir.pop();
iFlav=-1000; jFlav=-1000;
testp->Divide(xsec);
testp->Write("",TObject::kOverwrite);
pdfDir.pop();
}
fout->Write();
fout->Close();
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:16 PM (19 h, 39 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022918
Default Alt Text
(32 KB)

Event Timeline