Page MenuHomeHEPForge
Files F11222145

No OneTemporary
Actions

View Options

Index: example/trunk/example.cxx
===================================================================
--- example/trunk/example.cxx (revision 0)
+++ example/trunk/example.cxx (revision 1273)
@@ -0,0 +1,533 @@
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "appl_grid/appl_grid.h"
+#include "appl_grid/appl_pdf.h"
+
+
+#include "appl_grid/appl_timer.h"
+
+#include "TH1D.h"
+#include "TLatex.h"
+#include "TFile.h"
+#include "TH2D.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TStyle.h"
+#include "TColor.h"
+#include "TPaveStats.h"
+
+
+#include "appl_grid/lumi_pdf.h"
+
+
+#include "LHAPDF/PDFInterface.h"
+
+// lhapdf routines
+#include "LHAPDF/LHAPDF.h"
+extern "C" void evolvepdf_(const double& x, const double& Q, double* xf);
+extern "C" double alphaspdf_(const double& Q);
+
+
+
+extern "C" void mypdf(const double& x, const double& Q, double* xf) {
+
+ static double A = 208;
+ static double Z = 82;
+
+ evolvepdf_( x, Q, xf );
+
+ double u = Z*xf[7]/A + (A-Z)*xf[6]/A;
+ double d = Z*xf[6]/A + (A-Z)*xf[7]/A;
+
+ xf[7] = u;
+ xf[6] = d;
+
+ return;
+}
+
+
+
+void setpalette()
+{
+ const Int_t NRGBs = 7;
+ const Int_t NCont = 98;
+
+ Double_t stops[NRGBs] = { 0.00, 0.20, 0.40, 0.60, 0.70, 0.90, 1.00 };
+
+ Double_t red[NRGBs] = { 0.00, 0.00, 0.00, 0.60, 0.90, 1.00, 0.50 };
+ Double_t green[NRGBs] = { 0.00, 0.00, 0.40, 0.95, 1.00, 0.20, 0.00 };
+ Double_t blue[NRGBs] = { 0.00, 1.00, 1.00, 0.10, 0.00, 0.00, 0.00 };
+
+ TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+ gStyle->SetNumberContours(NCont);
+}
+
+
+
+
+
+TH2D* smooth(TH2D* h) {
+
+ int Nnbins = 10;
+
+ TH2D* hi = new TH2D("h2","",
+ h->GetXaxis()->GetNbins()*Nnbins,
+ h->GetXaxis()->GetBinLowEdge(1),
+ h->GetXaxis()->GetBinLowEdge(h->GetXaxis()->GetNbins()+1),
+
+ h->GetYaxis()->GetNbins()*Nnbins,
+ h->GetYaxis()->GetBinLowEdge(1),
+ h->GetYaxis()->GetBinLowEdge(h->GetYaxis()->GetNbins()+1) );
+
+ hi->GetXaxis()->SetTitle( h->GetXaxis()->GetTitle() );
+ hi->GetYaxis()->SetTitle( h->GetYaxis()->GetTitle() );
+ hi->GetYaxis()->SetTitleOffset(1.1);
+
+
+ for ( int i=1 ; i<=h->GetXaxis()->GetNbins() ; i++ ) {
+ for ( int j=1 ; j<=h->GetYaxis()->GetNbins() ; j++ ) {
+
+ double z1 = h->GetBinContent(i,j);
+ double z2 = h->GetBinContent(i+1,j);
+ double z3 = h->GetBinContent(i,j+1);
+ double z4 = h->GetBinContent(i+1,j+1);
+
+
+ if ( z2==0 ) z2 = z1;
+ if ( z3==0 ) z3 = z1;
+ if ( z4==0 ) z4 = z3;
+
+ for ( int ii=0 ; ii<Nnbins ; ii++ ) {
+
+ double fx = 1.0*ii/Nnbins;
+
+ for ( int jj=0 ; jj<Nnbins ; jj++ ) {
+
+ double fy = 1.0*jj/Nnbins;
+
+ double z = ( (z4-z3-z2+z1)*fx + z3 - z1 )*fy + (z2-z1)*fx + z1;
+
+ // std::cout << (i-1)*Nnbins+ii+1 << " " << (j-1)*Nnbins+jj+1 << "\t z : " << z << "\tz1-4: " << z1 << " " << z2 << " " << z3 << " " << z4 << std::endl;
+
+ hi->SetBinContent( (i-1)*Nnbins+ii+1, (j-1)*Nnbins+jj+1, z);
+ // hi->SetBinContent( (i-1)*Nnbins+ii+1, (j-1)*Nnbins+jj+1, h->GetBinContent(i,j) );
+
+ }
+
+ }
+
+ }
+
+ }
+
+ return hi;
+
+}
+
+
+
+
+int main(int argc, char** argv) {
+
+ /// which example to run
+ int example = 0;
+
+ /// use a default atlas inclusive grid
+ std::string gridname = "atlas-incljets-arxiv-1009.5908v2/r04/atlas-incljets-eta1.root";
+
+ if ( argc>1 ) {
+ for ( int i=1 ; i<argc ; i++ ) {
+ if ( std::string(argv[i])=="-e" || std::string(argv[i])=="--example" ) {
+ if ( i++<argc ) example = std::atoi(argv[1]);
+ else std::cerr << "need to give an example" << std::endl;
+ }
+ else gridname = argv[1];
+ }
+ }
+
+ gStyle->SetOptStat(0);
+
+ setpalette();
+
+
+
+ std::cout << "reading grid " << gridname << std::endl;
+
+
+ // initialise lhapdf
+
+ std::cout << "setting up lhapdf" << std::endl;
+
+ // const std::string _pdfname = "cteq6mE.LHgrid";
+ // const std::string _pdfname = "CT10.LHgrid";
+ const std::string _pdfname = "CT10nlo.LHgrid";
+ int iset = 1;
+
+ LHAPDF::initPDFSet( _pdfname, iset );
+ // initpdfset_(_pdfname.c_str());
+ // initpdf_( iset );
+
+
+
+ std::cout << "starting ... " << std::endl;
+
+ std::cout << "reading grid ... " << std::endl;
+
+ appl::grid g(gridname);
+
+ std::cout << "read" << std::endl;
+
+ /// ----------------- example 1 --------------------
+
+ if ( example==1 ) {
+
+ /// print the grid documentation
+ std::cout << g.getDocumentation() << std::endl;
+
+ // g.print();
+
+ // do the convolution into a vector
+ std::cout << "example 1: standalone convolution" << std::endl;
+
+ struct timeval atimer = appl_timer_start();
+
+ // std::vector<double> xsec = g.vconvolute( evolvepdf_, alphaspdf_ );
+
+ TH1D* hxsec = g.convolute( evolvepdf_, alphaspdf_ );
+
+
+ // std::vector<double> _xsec_duff = g.vconvolute( pdf, pdf );
+
+ // return 0;
+
+ double atime = appl_timer_stop(atimer);
+
+ std::cout << "time : " << atime << " ms" << std::endl;
+
+ atimer = appl_timer_start();
+
+ std::vector<double> xsec = g.vconvolute( evolvepdf_, alphaspdf_ );
+
+ atime = appl_timer_stop(atimer);
+
+ std::cout << "time : " << atime << " ms" << std::endl;
+
+ for ( int i=0 ; i<xsec.size() ; i++ ) std::cout << "xsec[" << i << "]\t= " << xsec[i] << std::endl;
+
+ TH1D* h = g.getReference();
+
+ for ( int i=0 ; i<xsec.size() ; i++ ) std::cout << i << " " << h->GetBinContent(i) << std::endl;
+
+
+
+ hxsec->SetTitle("");
+ // if ( i==0 )
+ hxsec->DrawCopy();
+ hxsec->DrawCopy("same");
+
+ // gPad->SetLogy(true);
+ gPad->Print("example-1.pdf");
+
+
+ TFile f("fout.root", "recreate");
+ hxsec->Write();
+ f.Write();
+ f.Close();
+
+ // }
+ }
+
+
+ /// ----------------- example 2 (corrections) --------------------
+
+ if ( example==2 ) {
+
+ std::cout << "example 2: multiplicative corrections" << std::endl;
+
+ std::cout << g.getApplyCorrections() << std::endl;
+
+ g.setApplyCorrections(false);
+ TH1D* huncor = g.convolute( evolvepdf_, alphaspdf_ );
+
+ g.setApplyCorrections(true);
+ TH1D* hcor = g.convolute( evolvepdf_, alphaspdf_ );
+
+ huncor->SetLineStyle(2);
+ huncor->SetTitle("");
+
+ huncor->DrawCopy();
+ hcor->DrawCopy("same");
+
+ gPad->Print("example-2.pdf");
+
+ }
+
+ /// ----------------- example 3 (cms scaling) --------------------
+
+ if ( example==3 ) {
+
+ std::cout << "example 3: centre-of-mass energy scaling " << std::endl;
+
+ TH1D* h7 = g.convolute( evolvepdf_, alphaspdf_ );
+
+ double Escale = 7/2.46;
+
+ TH1D* h246 = g.convolute( Escale, evolvepdf_, alphaspdf_ );
+
+ std::cout << "cms scale " << g.getCMSScale() << std::endl;
+
+ h7->SetTitle("");
+ h7->SetMinimum(0.001);
+
+ gStyle->SetOptStat(0);
+
+ h7->DrawCopy();
+ h246->SetLineColor(kRed);
+ h246->DrawCopy("same");
+
+ gPad->SetLogy(true);
+ gPad->Print("example-3.pdf");
+
+ }
+
+
+
+
+
+ /// ------------------ example 4 (scale variation) -----------------
+
+ if ( example==4 ) {
+
+ std::cout << "example 4: scale variation " << std::endl;
+
+ int nloops = 1; // loop order
+
+ int Nbins = 20;
+ TH2D* hscale = new TH2D("scale", "", Nbins, 0.45, 2.05, Nbins, 0.45, 2.05 );
+
+ struct timeval atimer = appl_timer_start();
+
+ for ( int i=1 ; i<=hscale->GetNbinsX() ; i++ ) {
+ double fscale = hscale->GetXaxis()->GetBinCenter(i);
+
+ std::cout << i << "\tfscale " << fscale << std::endl;
+
+ for ( int j=1 ; j<=hscale->GetNbinsY() ; j++ ) {
+ double rscale = hscale->GetYaxis()->GetBinCenter(j);
+
+ std::vector<double> xs = g.vconvolute( evolvepdf_, alphaspdf_, nloops, rscale, fscale );
+
+ double total = 0;
+ for ( int k=xs.size() ; k-- ; ) total += xs[k];
+
+ hscale->Fill( fscale, rscale, total );
+
+ }
+ }
+
+ double atime = appl_timer_stop(atimer);
+
+ std::cout << "time " << atime << " ms" << std::endl;
+
+
+ /// plot ....
+
+
+ gStyle->SetPalette(52);
+ gStyle->SetPadRightMargin(0.17);
+
+
+ TCanvas* tc = new TCanvas( "","", 650, 600 );
+
+ hscale->GetXaxis()->SetTitle("factorisation scale factor");
+ hscale->GetYaxis()->SetTitle("renormalisation scale factor");
+ hscale->GetYaxis()->SetTitleOffset(1.1);
+
+ TH2D* _hscale = smooth(hscale);
+
+ _hscale->DrawCopy("colz");
+
+ gPad->Print("example-4.png");
+
+
+ TFile fs("scale.root","recreate");
+ fs.cd();
+ hscale->Write();
+ fs.Write();
+ fs.Close();
+
+ }
+
+
+
+ /// ------------------ example 5 (subprocesses) -----------------
+
+ if ( example==5 ) {
+
+ std::cout << "example 4: subprocess contributions " << std::endl;
+
+ std::vector<double> xsec = g.vconvolute( evolvepdf_, alphaspdf_);
+
+ std::vector<TH1D*> xsec_sub(g.subProcesses(0));
+
+ for ( int i=0 ; i<g.subProcesses(0) ; i++ ) {
+ xsec_sub[i] = g.convolute_subproc( i, evolvepdf_, alphaspdf_);
+
+ for ( int j=0 ; j<xsec_sub[i]->GetNbinsX() ; j++ ) {
+ xsec_sub[i]->SetBinContent( j+1, xsec_sub[i]->GetBinContent(j+1)/xsec[j] );
+ }
+ }
+
+
+ /// plot ...
+
+ int col[7] = { kRed, kBlue, kMagenta, kCyan, kGreen, kBlack, kViolet };
+
+ xsec_sub[0]->SetMaximum(0.6);
+ xsec_sub[0]->SetMinimum(0);
+ xsec_sub[0]->GetXaxis()->SetTitle("P_{T} [GeV]");
+ xsec_sub[0]->GetYaxis()->SetTitle("subprocess fraction");
+ xsec_sub[0]->SetTitle("");
+
+
+ xsec_sub[0]->DrawCopy();
+
+
+ for ( int i=0 ; i<g.subProcesses(0) ; i++ ) {
+ xsec_sub[i]->SetLineColor(col[i]);
+ xsec_sub[i]->DrawCopy("same");
+
+ TLatex* tl = new TLatex(0,0,"");
+ tl->SetTextColor(col[i]);
+ tl->SetTextSize(0.04);
+
+ char label[64];
+ sprintf(label, "subprocess %d", i);
+
+ tl->DrawLatex(430, 0.55-0.035*i, label );
+ }
+
+ gPad->SetLogy(false);
+
+ gPad->Print("example-5.pdf");
+
+ }
+
+
+
+
+ /// -------------------- example 6 (pdf sets) -----------------------
+
+ if ( example==6 ) {
+
+ // const std::string _pdfname = "cteq6mE.LHgrid";
+ const std::string _pdfname = "CT10.LHgrid";
+ int iset = 0;
+
+ LHAPDF::initPDFSet( _pdfname, iset );
+ std::vector<double> _xs = g.vconvolute( evolvepdf_, alphaspdf_ );
+
+ int Nset = LHAPDF::numberPDF();
+
+ std::vector<TH1D*> xs(Nset);
+
+ for ( int i=0 ; i<Nset ; i++ ) {
+ LHAPDF::initPDFSet( _pdfname, i );
+ xs[i] = g.convolute( evolvepdf_, alphaspdf_ );
+
+ for ( int j=_xs.size() ; j-- ; ) xs[i]->SetBinContent( j+1, xs[i]->GetBinContent(j+1)/_xs[j]);
+ }
+
+ xs[0]->SetMaximum(1.1);
+ xs[0]->SetMinimum(0.9);
+ gPad->SetLogy(false);
+ gPad->SetLogx(true);
+
+ xs[0]->DrawCopy("l");
+ for ( int i=0 ; i<Nset ; i++ ) xs[i]->DrawCopy("lsame");
+
+ gPad->Print("example-6.pdf");
+ gPad->SetLogx(false);
+ }
+
+
+ /// ------------------- example 7 (timers) ---------------------
+
+ if ( example==7 ) {
+
+ gStyle->SetOptStat(2220);
+
+ std::cout << "example 7 : timings" << std::endl;
+
+ TH1D* htrim = new TH1D("htrim", "", 100, 38, 45);
+ // TH1D* htrim2 = new TH1D("htrim2", "", 100, 37, 45);
+ TH1D* huntrim = new TH1D("huntrim", "", 100, 38, 45);
+
+
+ g.trim();
+
+ std::cout << "grid size : " << g.size() << std::endl;
+
+ for ( int i=400 ; i-- ; ) {
+ struct timeval atimer = appl_timer_start();
+ std::vector<double> xsec = g.vconvolute( evolvepdf_, alphaspdf_ );
+ double atime = appl_timer_stop( atimer );
+ htrim->Fill(atime);
+ }
+
+ g.untrim();
+
+ std::cout << "grid size : " << g.size() << std::endl;
+
+ for ( int i=400 ; i-- ; ) {
+ struct timeval atimer = appl_timer_start();
+ std::vector<double> xsec = g.vconvolute( evolvepdf_, alphaspdf_ );
+ double atime = appl_timer_stop( atimer );
+ huntrim->Fill(atime);
+ }
+
+
+ if ( huntrim->GetMaximum()<htrim->GetMaximum() ) huntrim->SetMaximum(1.1*htrim->GetMaximum());
+ else htrim->SetMaximum(1.1*huntrim->GetMaximum());
+ huntrim->SetLineStyle(2);
+ huntrim->SetLineColor(kRed);
+ htrim->GetXaxis()->SetTitle("time [ms]");
+ huntrim->GetXaxis()->SetTitle("time [ms]");
+
+ huntrim->Draw();
+ htrim->Draw("same");
+ // htrim2->DrawCopy("same");
+ huntrim->Draw("same");
+ htrim->Draw("sames+");
+ huntrim->Draw("sames+");
+
+ gPad->Update();
+
+ TPaveStats* tp1 = (TPaveStats*)htrim->GetListOfFunctions()->FindObject("stats");
+ TPaveStats* tp2 = (TPaveStats*)huntrim->GetListOfFunctions()->FindObject("stats");
+
+ std::cout << tp1 << " " << tp2 << std::endl;
+
+ tp2->SetTextColor(kRed);
+
+
+ tp1->SetY1NDC(0.93-0.15);
+ tp1->SetY2NDC(0.93);
+
+ tp2->SetY1NDC(0.93-0.16-0.15);
+ tp2->SetY2NDC(0.93-0.16);
+
+ tp1->Draw();
+ tp2->Draw();
+
+ gPad->Update();
+
+ gPad->Print("example-7.pdf");
+
+ }
+
+ return 0;
+}
Index: example/trunk/version.dat
===================================================================
--- example/trunk/version.dat (revision 1272)
+++ example/trunk/version.dat (revision 1273)
@@ -1 +1 @@
-0.0.2
+0.0.3
Index: example/trunk/makefile
===================================================================
--- example/trunk/makefile (revision 1272)
+++ example/trunk/makefile (revision 1273)
@@ -1,86 +1,89 @@
CXX = g++
F77 = gfortran
FFLAGS += -O3 -fPIC
CXXFLAGS += -O3 -fPIC
# root
ROOTINCS = $(shell root-config --cflags)
ROOTLIBS = $(shell root-config --glibs)
ROOTARCH = $(findstring -m64, $(ROOTINCS) )
#LHAPDF
LHAPDFINCS = -I$(shell lhapdf-config --prefix)/include
LHAPDFDIR = $(shell lhapdf-config --prefix)/lib
LHAPDFLIBS = -L$(LHAPDFDIR) -lLHAPDF
# applgrid
APPLCXXFLAGS = $(shell applgrid-config --cxxflags)
APPLCLIBS = $(shell applgrid-config --ldflags)
APPLFLIBS = $(shell applgrid-config --ldfflags)
# hoppet
HOPPETLIBS = $(shell hoppet-config --libs)
# get the fotran runtime library for linking fortran
FRTLLIB = $(shell gfortran $(CXXFLAGS) -print-file-name=libgfortran.a)
FRTLIB = -L$(subst /libgfortran.a, ,$(FRTLLIB) ) -lgfortran
# now set up the compile and link flags and libs
CXXFLAGS += $(ROOTARCH) $(ROOTINCS) $(APPLCXXFLAGS) $(LHAPDFINCS)
LDFLAGS += $(ROOTARCH)
FFLAGS += $(ROOTARCH)
CLIBS += $(ROOTLIBS) $(LHAPDFLIBS) $(HOPPETLIBS) $(APPLCLIBS)
FLIBS += $(ROOTLIBS) $(LHAPDFLIBS) $(HOPPETLIBS) $(APPLFLIBS) $(APPLCLIBS) $(FRTLIB)
install : all
-all : stand fnlo fstand getdata
+all : stand fnlo fstand example getdata
stand: stand.o
$(CXX) $(LDFLAGS) -o $@ $< $(CLIBS)
fstand: fstand.o fmain.o
$(CXX) $(LDFLAGS) -o $@ fstand.o fmain.o $(FLIBS)
+example: example.o
+ $(CXX) $(LDFLAGS) -o $@ $< $(CLIBS)
+
% : %.o
$(CXX) $(LDFLAGS) -o $@ $< $(CLIBS)
fnlo: fnmain.o
$(CXX) $(LDFLAGS) -o $@ $< $(CLIBS)
.SUFFIXES : .cxx .o .f .c
.f.o :
$(F77) $(FFLAGS) -c $<
.cxx.o:
$(CXX) $(CXXFLAGS) -c $<
#########################################################
# just type make getdata and it will download all
# the required data for these examples
#########################################################
.PRECIOUS : fnlodata appldata
getdata : getfnlodata getappldata
getfnlodata: fnl0004.tab
fnl0004.tab:
@ echo "\n************************\ndownloading fastnlo data\n************************\n"
(wget http://fastnlo.hepforge.org/code/src/fnl0004.tab.gz ; gunzip fnl0004.tab.gz > fnl0004.tab )
getappldata: atlas-incljets-arxiv-1009.5908v2
atlas-incljets-arxiv-1009.5908v2:
@ echo "\n*************************************\n downloading atlas inclusive jet data\n***********************************\n"
(wget http://www.hepforge.org/archive/applgrid/atlas-incljets-arxiv-1009.5908v2.tgz ; tar -xzf atlas-incljets-arxiv-1009.5908v2.tgz )
clean:
rm -rf ./.libs ./.obj *.lo *.o *.la stand fnlo fstand *~

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:29 AM (11 h, 22 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111435
Default Alt Text
(15 KB)

Event Timeline