diff --git a/Cards/Config/Core.py b/Cards/Config/Core.py index 76b0a9c..8396363 100644 --- a/Cards/Config/Core.py +++ b/Cards/Config/Core.py @@ -1,165 +1,168 @@ #!/usr/bin/env python class PrintHelper(object): _indent = 0 _indent_size = 4 def indent(self): self._indent += self._indent_size def unindent(self): self._indent -= self._indent_size def indentation(self): return ' '*self._indent class Parameters(dict): '''A raw list of steering parameters''' __getattr__ = dict.get __setattr__ = dict.__setitem__ __delattr__ = dict.__delitem__ def __init__(self, *args, **kwargs): self.update(*args, **kwargs) super(Parameters, self).__init__(*args, **kwargs) def __deepcopy__(self, memo): from copy import deepcopy return Parameters([(deepcopy(k, memo), deepcopy(v, memo)) for k, v in self.items()]) def dump(self, printer=PrintHelper()): out = self.__class__.__name__+'(\n' for k, v in self.items(): printer.indent() out += ('%s%s = ' % (printer.indentation(), k)) if v.__class__.__name__ not in ['Parameters', 'Module']: out += v.__repr__() else: out += v.dump(printer) out += ',\n' printer.unindent() out += printer.indentation()+')' return out def __repr__(self): return self.dump() def clone(self, *args, **kwargs): from copy import deepcopy out = deepcopy(self) for k in kwargs: out[k] = kwargs.get(k) return type(self)(out) def load(self, mod): mod = mod.replace('/', '.') module = __import__(mod) self.extend(sys.modules[mod]) class Module(Parameters): '''A named parameters set to steer a generic module''' def __init__(self, name, *args, **kwargs): super(Module, self).__init__(*args, **kwargs) self.mod_name = name def __len__(self): return dict.__len__(self)-1 # discard the name key def dump(self, printer=PrintHelper()): out = self.__class__.__name__+'('+self.mod_name.__repr__()+'\n' mod_repr = self.clone('') mod_repr.pop('mod_name', None) for k, v in mod_repr.items(): printer.indent() out += ('%s%s = ' % (printer.indentation(), k)) if v.__class__.__name__ not in ['Parameters', 'Module']: out += v.__repr__() else: out += v.dump(printer) out += ',\n' printer.unindent() out += printer.indentation()+')' return out def __repr__(self): return self.dump() def clone(self, name, **kwargs): out = Parameters(self).clone(**kwargs) out.mod_name = name return out class Logging: '''Logging verbosity''' Nothing = 0 Error = 1 Warning = 2 Information = 3 Debug = 4 DebugInsideLoop = 5 class StructureFunctions: '''Types of structure functions supported''' Electron = Parameters(id=1) ElasticProton = Parameters(id=2) SuriYennie = Parameters(id=11) SzczurekUleshchenko = Parameters(id=12) BlockDurandHa = Parameters(id=13) FioreBrasse = Parameters(id=101) ChristyBosted = Parameters(id=102) CLAS = Parameters(id=103) ALLM91 = Parameters(id=201) ALLM97 = Parameters(id=202) GD07p = Parameters(id=203) GD11p = Parameters(id=204) - MSTWgrid = Parameters(id=205) + MSTWgrid = Parameters( + id = 205, + gridPath = 'External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat', + ) LUXlike = Parameters( id = 301, #Q2cut = 10., #W2limits = (4.,1.), continuumSF = GD11p, resonancesSF = ChristyBosted, ) GenericLHAPDF = Parameters( id = 401, pdfSet = 'LUXqed17_plus_PDF4LHC15_nnlo_100', numFlavours = 5, ) class ProcessMode: '''Types of processes supported''' ElectronProton = 0 ElasticElastic = 1 ElasticInelastic = 2 InelasticElastic = 3 InelasticInelastic = 4 ProtonElectron = 5 ElectronElectron = 6 if __name__ == '__main__': import unittest class TestTypes(unittest.TestCase): def testModules(self): mod = Module('empty') self.assertEqual(len(mod), 0) mod.param1 = 'foo' self.assertEqual(len(mod), 1) # playing with modules clones mod_copy = mod.clone('notEmpty', param1 = 'boo', param2 = 'bar') self.assertEqual(mod.param1, 'foo') self.assertEqual(mod_copy.param1, 'boo') self.assertEqual(mod_copy.param2, 'bar') self.assertEqual(mod.param1+mod_copy.param2, 'foobar') def testParameters(self): params = Parameters( first = 'foo', second = 'bar', third = 42, fourth = (1, 2), ) params_copy = params.clone( second = 'bak', ) self.assertEqual(len(params), 4) self.assertEqual(params.first, params['first']) self.assertEqual(params['second'], 'bar') self.assertTrue(int(params.third) == params.third) self.assertEqual(len(params.fourth), 2) self.assertEqual(params.second, 'bar') # playing with parameters clones self.assertEqual(params_copy.second, 'bak') # check that the clone does not change value if the origin does # (i.e. we indeed have a deep copy and not a shallow one...) params.third = 43 self.assertEqual(params.third, 43) self.assertEqual(params_copy.third, 42) unittest.main() diff --git a/CepGen/StructureFunctions/MSTWGrid.cpp b/CepGen/StructureFunctions/MSTWGrid.cpp index 1d7292a..e73648f 100644 --- a/CepGen/StructureFunctions/MSTWGrid.cpp +++ b/CepGen/StructureFunctions/MSTWGrid.cpp @@ -1,119 +1,118 @@ #include "CepGen/StructureFunctions/MSTWGrid.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include <fstream> -#include <atomic> #include <gsl/gsl_errno.h> #include <gsl/gsl_math.h> namespace MSTW { const unsigned int Grid::good_magic = 0x5754534d; // MSTW in ASCII Grid& Grid::get( const char* filename ) { Parameterisation p; p.grid_path = filename; static Grid instance( p ); return instance; } Grid::Grid( const Parameterisation& param ) : CepGen::StructureFunctions( CepGen::SF::Type::MSTWgrid ), CepGen::GridHandler<2>( CepGen::GridType::kLogarithmic ), params( param ) { std::set<double> q2_vals, xbj_vals; { // file readout part std::ifstream file( params.grid_path, std::ios::binary | std::ios::in ); if ( !file.is_open() ) throw CG_FATAL( "Grid" ) << "Impossible to load grid file \"" << params.grid_path << "\"!"; file.read( reinterpret_cast<char*>( &header_ ), sizeof( header_t ) ); // first checks on the file header if ( header_.magic != good_magic ) throw CG_FATAL( "Grid" ) << "Wrong magic number retrieved: " << header_.magic << ", expecting " << good_magic << "."; if ( header_.nucleon != header_t::proton ) throw CG_FATAL( "Grid" ) << "Only proton structure function grids can be retrieved for this purpose!"; // retrieve all points and evaluate grid boundaries sfval_t val; while ( file.read( reinterpret_cast<char*>( &val ), sizeof( sfval_t ) ) ) { q2_vals.insert( log10( val.x ) ); xbj_vals.insert( log10( val.y ) ); #ifndef GOOD_GSL x_vals_.emplace_back( val.x ); y_vals_.emplace_back( val.y ); #endif values_raw_.emplace_back( val ); } file.close(); } if ( q2_vals.size() < 2 || xbj_vals.size() < 2 ) throw CG_FATAL( "Grid" ) << "Invalid grid retrieved!"; initGSL( q2_vals, xbj_vals ); CG_INFO( "Grid" ) << "MSTW@" << header_.order << " grid evaluator built " << "for " << header_.nucleon << " structure functions (" << header_.cl << ")\n\t" << " Q² in range [" << pow( 10., *q2_vals.begin() ) << ":" << pow( 10., *q2_vals.rbegin() ) << "]\n\t" << "xBj in range [" << pow( 10., *xbj_vals.begin() ) << ":" << pow( 10., *xbj_vals.rbegin() ) << "]."; } Grid& Grid::operator()( double q2, double xbj ) { const sfval_t val = CepGen::GridHandler<2>::eval( q2, xbj ); F2 = val.value[0]; FL = val.value[1]; return *this; } std::ostream& operator<<( std::ostream& os, const sfval_t& val ) { return os << Form( "Q² = %.5e GeV²\txbj = %.4f\tF₂ = % .6e\tFL = % .6e", val.x, val.y, val.value[0], val.value[1] ); } std::ostream& operator<<( std::ostream& os, const Grid::header_t::order_t& order ) { switch ( order ) { case Grid::header_t::lo: return os << "LO"; case Grid::header_t::nlo: return os << "nLO"; case Grid::header_t::nnlo: return os << "nnLO"; } return os; } std::ostream& operator<<( std::ostream& os, const Grid::header_t::cl_t& cl ) { switch ( cl ) { case Grid::header_t::cl68: return os << "68% C.L."; case Grid::header_t::cl95: return os << "95% C.L."; } return os; } std::ostream& operator<<( std::ostream& os, const Grid::header_t::nucleon_t& nucl ) { switch ( nucl ) { case Grid::header_t::proton: return os << "proton"; case Grid::header_t::neutron: return os << "neutron"; } return os; } } diff --git a/test/Canvas.h b/test/Canvas.h index ecf674a..13e912a 100644 --- a/test/Canvas.h +++ b/test/Canvas.h @@ -1,320 +1,318 @@ #ifndef Canvas_h #define Canvas_h #include "TCanvas.h" #include "TLegend.h" #include "TPaveText.h" #include "TObjArray.h" #include "TObjString.h" #include "TH1.h" #include "TGraphErrors.h" #include "TStyle.h" #include <string.h> #define font_type(x) 130+x namespace CepGen { class PaveText : public TPaveText { public: inline PaveText( const float& x1, const float& y1, const float& x2, const float& y2, const char* text="" ) : TPaveText( x1, y1, x2, y2, "NDC" ) { TPaveText::SetTextAlign( 13 ); if ( strcmp( text, "" )!=0 ) { TString txt = text; if ( txt.Contains( "\\" ) ) { TObjArray* tok = txt.Tokenize( "\\" ); for ( int i=0; i<tok->GetEntries(); i++ ) { TPaveText::AddText( dynamic_cast<TObjString*>( tok->At( i ) )->String() ); } } else TPaveText::AddText( text ); } TPaveText::SetFillColor( 0 ); TPaveText::SetFillStyle( 0 ); TPaveText::SetLineColor( 0 ); TPaveText::SetLineWidth( 0 ); TPaveText::SetShadowColor( 0 ); TPaveText::SetTextFont( font_type( 2 ) ); TPaveText::SetTextSize( 0.058 ); } }; class Canvas : public TCanvas { public: inline Canvas( const char* name, const char* title="", bool ratio=false ) : //TCanvas( name, "", 450, 450 ), TCanvas( name, "", 600, 600 ), fTitle( title ), fTopLabel( 0 ), fLeg( 0 ), fLegX1( 0.5 ), fLegY1( 0.75 ), fRatio( ratio ) { gStyle->SetOptStat( 0 ); Build(); } inline ~Canvas() { if ( fLeg ) delete fLeg; if ( fTopLabel ) delete fTopLabel; } inline void SetSize( const float& size=600 ) { TCanvas::SetCanvasSize( size, 600 ); } inline void Prettify( TH1* obj ) { TAxis* x = dynamic_cast<TAxis*>( obj->GetXaxis() ), *y = dynamic_cast<TAxis*>( obj->GetYaxis() ), *z = dynamic_cast<TAxis*>( obj->GetZaxis() ); x->SetLabelFont( font_type( 3 ) ); x->SetLabelSize( 20 ); x->SetTitleFont( font_type( 3 ) ); x->SetTitleSize( 29 ); y->SetLabelFont( font_type( 3 ) ); y->SetLabelSize( 20 ); y->SetTitleFont( font_type( 3 ) ); y->SetTitleSize( 29 ); z->SetLabelFont( font_type( 3 ) ); z->SetLabelSize( 16 ); z->SetTitleFont( font_type( 3 ) ); z->SetTitleSize( 29 ); if ( fRatio ) { x->SetTitleOffset( 3. ); x->SetLabelOffset( 0.02 ); } y->SetTitleOffset( 1.3 ); x->SetTickLength( 0.03 ); y->SetTickLength( 0.03 ); // axis titles TString ttle = obj->GetTitle(); if ( ttle.Contains( "\\" ) ) { TObjArray* tok = ttle.Tokenize( "\\" ); TString x_title = "", y_title = "", unit = "", form_spec = "", distrib = ""; if ( tok->GetEntries()>0 ) x_title = dynamic_cast<TObjString*>( tok->At( 0 ) )->String(); if ( tok->GetEntries()>1 ) y_title = dynamic_cast<TObjString*>( tok->At( 1 ) )->String(); if ( tok->GetEntries()>2 ) { unit = ( ( TObjString* )tok->At( 2 ) )->String(); if ( unit.Contains( "?" ) ) { // extract format specifier TObjArray* tok2 = unit.Tokenize( "?" ); if ( tok2->GetEntries()>1 ) { unit = dynamic_cast<TObjString*>( tok2->At( 0 ) )->String(); form_spec = dynamic_cast<TObjString*>( tok2->At( 1 ) )->String(); } else { unit = ""; form_spec = dynamic_cast<TObjString*>( tok2->At( 0 ) )->String(); } } } if ( tok->GetEntries()>3 ) { distrib = ( ( TObjString* )tok->At( 3 ) )->String(); } if ( !unit.IsNull() or !form_spec.IsNull() ) { if ( !unit.IsNull() ) x_title = Form( "%s (%s)", x_title.Data(), unit.Data() ); if ( !distrib.IsNull() ) { if ( !form_spec.IsNull() ) { TString format = Form( "%%s (%s / %%%s %%s)", distrib.Data(), form_spec.Data() ); y_title = Form( format.Data(), y_title.Data(), GetBinning( obj ), unit.Data() ); } else y_title = Form( "%s (%s / %d %s)", y_title.Data(), distrib.Data(), static_cast<unsigned int>( GetBinning( obj ) ), unit.Data() ); } - else { + else { if ( !form_spec.IsNull() ) { TString format = Form( "%%s / %%%s %%s", form_spec.Data() ); y_title = Form( format.Data(), y_title.Data(), GetBinning( obj ), unit.Data() ); } else y_title = Form( "%s / %d %s", y_title.Data(), static_cast<unsigned int>( GetBinning( obj ) ), unit.Data() ); } } obj->GetXaxis()->SetTitle( x_title ); obj->GetYaxis()->SetTitle( y_title ); obj->SetTitle( "" ); } //else obj->GetXaxis()->SetTitle(ttle); } inline void DrawDiagonal(const TH1* obj) { TLine l; l.SetLineWidth( 2 ); l.SetLineColor( kGray ); l.SetLineStyle( 2 ); l.DrawLine( obj->GetXaxis()->GetXmin(), obj->GetYaxis()->GetXmin(), obj->GetXaxis()->GetXmax(), obj->GetYaxis()->GetXmax() ); } inline void RatioPlot(TH1* obj1, const TH1* obj2, const TH1* obj3, float ymin=-999., float ymax=-999.) { if (!fRatio) return; TH1* ratio1 = (TH1*)obj2->Clone(), *ratio2 = (TH1*)obj3->Clone(); //ratio1->Sumw2(); ratio2->Sumw2(); ratio1->Divide(obj1); ratio2->Divide(obj1); TCanvas::cd(2); ratio1->Draw("p"); ratio2->Draw("p same"); obj1->GetXaxis()->SetTitle(""); if ( ymin!=ymax ) { ratio1->GetYaxis()->SetRangeUser(ymin, ymax); } Prettify(ratio1); TCanvas::cd(); } inline void RatioPlot( TH1* obj1, const TH1* obj2=0, float ymin=-999., float ymax=-999. ) { if ( !fRatio ) return; TH1* ratio; if ( obj2 ) { ratio = dynamic_cast<TH1*>( obj2->Clone() ); ratio->Divide( obj1 ); } else { ratio = dynamic_cast<TH1*>( obj1->Clone() ); } TCanvas::cd( 2 ); ratio->Draw("p"); obj1->GetXaxis()->SetTitle(""); if ( ymin!=ymax ) { ratio->GetYaxis()->SetRangeUser( ymin, ymax ); } Prettify(ratio); ratio->GetYaxis()->SetTitle("Ratio"); TCanvas::cd(); } inline TGraphErrors* RatioPlot(TGraphErrors* obj1, const TGraphErrors* obj2, float ymin=-999., float ymax=-999.) { if (!fRatio) return 0; TGraphErrors* ratio = new TGraphErrors; ratio->SetTitle( obj1->GetTitle() ); unsigned int n = 0; float min_x = 9.e10, max_x = -9.e10; for ( int i=0; i<obj1->GetN(); i++ ) { const float x1 = obj1->GetX()[i]; for ( int j=0; j<obj2->GetN(); j++ ) { const float x2 = obj2->GetX()[j]; if (x2>max_x) max_x = x2; if (x2<min_x) min_x = x2; if ( fabs( x2-x1 )>1.e-3 ) continue; const float y1 = obj1->GetY()[i], y1_err = obj1->GetEY()[i], y2 = obj2->GetY()[j], y2_err = obj2->GetEY()[j]; const float y = (y2-y1)/y1, err_y = sqrt( pow( y1_err/y1, 2 )+pow( y2_err/y2, 2 )*y2/y1 ); ratio->SetPoint( n, x1, y ); ratio->SetPointError( n, 0., err_y ); n++; } } TCanvas::cd(2); ratio->Draw("ap"); ratio->GetXaxis()->SetRangeUser( obj1->GetXaxis()->GetXmin(), obj1->GetXaxis()->GetXmax() ); ratio->SetMarkerStyle( 20 ); if ( ymin!=ymax ) { ratio->GetYaxis()->SetRangeUser(ymin, ymax); } ratio->GetXaxis()->SetLimits(min_x, max_x); Prettify( ratio->GetHistogram() ); obj1->GetXaxis()->SetTitle(""); TLine l( min_x, 0., max_x, 0. ); l.Draw(); ratio->GetYaxis()->SetLabelSize( 14 ); TCanvas::cd(); return ratio; } inline void SetTopLabel(const char* lab="") { TCanvas::cd(); if (strcmp(lab, "")!=0) fTitle = lab; if (!fTopLabel) BuildTopLabel(); else fTopLabel->Clear(); fTopLabel->AddText(fTitle); //fTopLabel->Draw(); } inline void SetLegendX1(double x) { fLegX1 = x; } inline void SetLegendY1(double y) { fLegY1 = y; } inline void AddLegendEntry(const TObject* obj, const char* title, Option_t* option="lpf") { if (!fLeg) BuildLeg(); fLeg->AddEntry(obj, title, option); const unsigned int num_entries = fLeg->GetNRows(); if ( num_entries>3 ) { fLeg->SetY1( fLeg->GetY1()-( num_entries-3 )*0.015 ); } } inline void Save(const char* ext, const char* out_dir=".") { - if (strstr(ext, "pdf")==NULL) { - if (strstr(ext, "png")==NULL) { - if (strstr(ext, "root")==NULL) { - return; - } - } - } + if (strstr(ext, "pdf")==NULL) + if (strstr(ext, "eps")==NULL) + if (strstr(ext, "png")==NULL) + if (strstr(ext, "root")==NULL) + return; TCanvas::cd(); if ( fLeg ) fLeg->Draw(); if ( fTopLabel ) fTopLabel->Draw(); TCanvas::SaveAs(Form("%s/%s.%s", out_dir, TCanvas::GetName(), ext)); } inline TLegend* GetLegend() { return fLeg; } private: inline void Build() { TCanvas::SetLeftMargin(0.14); TCanvas::SetTopMargin(0.06); TCanvas::SetRightMargin(0.1); TCanvas::SetBottomMargin(0.12); TCanvas::SetTicks(1,1); SetTopLabel(); if (fRatio) DivideCanvas(); } inline void DivideCanvas() { TCanvas::Divide(1,2); TPad* p1 = (TPad*)TCanvas::GetPad(1), *p2 = (TPad*)TCanvas::GetPad(2); p1->SetPad(0., 0.3, 1., 1.); p2->SetPad(0., 0.0, 1., 0.3); p1->SetLeftMargin(TCanvas::GetLeftMargin()); p1->SetRightMargin(TCanvas::GetRightMargin()); p2->SetLeftMargin(TCanvas::GetLeftMargin()); p2->SetRightMargin(TCanvas::GetRightMargin()); p1->SetTopMargin(TCanvas::GetTopMargin()+0.025); p1->SetBottomMargin(0.02); p2->SetTopMargin(0.02); p2->SetBottomMargin(TCanvas::GetBottomMargin()+0.25); p1->SetTicks(1,1); p2->SetTicks(1,1); p2->SetGrid(0,1); TCanvas::cd(1); } inline void BuildTopLabel() { TCanvas::cd(); fTopLabel = new TPaveText(0.5, 0.95, 0.915, 0.96, "NB NDC"); fTopLabel->SetFillStyle(0); fTopLabel->SetFillColor(0); fTopLabel->SetLineColor(0); fTopLabel->SetLineStyle(0); fTopLabel->SetTextFont( font_type( 2 ) ); fTopLabel->SetTextSize(0.04); fTopLabel->SetTextAlign(kHAlignRight+kVAlignBottom); } inline void BuildLeg() { if ( fLeg ) return; if ( fRatio ) TCanvas::cd(1); fLeg = new TLegend(fLegX1, fLegY1, fLegX1+0.3, fLegY1+0.15); fLeg->SetLineColor(kWhite); fLeg->SetLineWidth(0); fLeg->SetFillStyle(0); fLeg->SetTextFont( font_type( 2 ) ); fLeg->SetTextSize(0.04); } inline float GetBinning(const TH1* h) { return (h->GetXaxis()->GetXmax()-h->GetXaxis()->GetXmin())/h->GetXaxis()->GetNbins(); } TString fTitle; TPaveText* fTopLabel; TLegend* fLeg; double fLegX1, fLegY1; bool fRatio; }; } #endif diff --git a/test/test_physics_processes.cpp b/test/test_physics_processes.cpp index 5db4939..fecd592 100644 --- a/test/test_physics_processes.cpp +++ b/test/test_physics_processes.cpp @@ -1,206 +1,205 @@ #include "CepGen/Generator.h" #include "CepGen/Parameters.h" #include "CepGen/Core/Timer.h" #include "CepGen/Processes/GamGamLL.h" #include "CepGen/Processes/PPtoFF.h" #include "CepGen/Processes/PPtoWW.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/PDG.h" #include "abort.h" #include <unordered_map> #include <assert.h> #include <string.h> using namespace std; int main( int argc, char* argv[] ) { typedef vector<pair<const char*,pair<double,double> > > KinematicsMap; typedef vector<pair<float, KinematicsMap> > ValuesAtCutMap; AbortHandler ctrl_c; // values defined at pt(single lepton)>15 GeV, |eta(single lepton)|<2.5, mX<1000 GeV // process -> { pt cut -> { kinematics -> ( sigma, delta(sigma) ) } } vector<pair<const char*,ValuesAtCutMap> > values_map = { //--- LPAIR values at sqrt(s) = 13 TeV /*{ "lpair", { { 3.0, { // pt cut { "elastic", { 2.0871703e1, 3.542e-2 } }, { "singlediss", { 1.5042536e1, 3.256e-2 } }, { "doublediss", { 1.38835e1, 4.03018e-2 } } } }, { 15.0, { // pt cut { "elastic", { 4.1994803e-1, 8.328e-4 } }, { "singlediss", { 4.8504819e-1, 1.171e-3 } }, { "doublediss", { 6.35650e-1, 1.93968e-3 } } } }, } },*/ //--- PPtoLL values { "pptoll", { { 3.0, { // pt cut { "elastic", { 2.0936541e1, 1.4096e-2 } }, { "singlediss_su", { 1.4844881e1, 2.0723e-2 } }, // SU, qt<50 { "doublediss_su", { 1.3580637e1, 2.2497e-2 } }, // SU, qt<50 } }, { 15.0, { // pt cut { "elastic", { 4.2515888e-1, 3.0351e-4 } }, { "singlediss_su", { 4.4903253e-1, 5.8970e-4 } }, // SU, qt<50 { "doublediss_su", { 5.1923819e-1, 9.6549e-4 } }, // SU, qt<50 /*{ "2_singlediss", { 4.6710287e-1, 6.4842e-4 } }, // SU, qt<500 { "3_doublediss", { 5.6316353e-1, 1.1829e-3 } }, // SU, qt<500*/ } }, } }, //--- PPtoWW values { "pptoww", { { 0.0, { // pt cut { "elastic", { 0.273, 0.01 } }, { "elastic", { 0.273, 0.01 } }, // FIXME { "singlediss_lux", { 0.409, 0.01 } }, { "doublediss_lux", { 1.090, 0.01 } }, { "singlediss_allm", { 0.318, 0.01 } }, { "doublediss_allm", { 0.701, 0.01 } } } } } }, }; const double num_sigma = 3.0; - if ( argc < 3 || strcmp( argv[2], "debug" ) != 0 ) { + if ( argc < 3 || strcmp( argv[2], "debug" ) != 0 ) CepGen::Logger::get().level = CepGen::Logger::Level::nothing; - } CepGen::Timer tmr; CepGen::Generator mg; if ( argc > 1 && strcmp( argv[1], "plain" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::Type::plain; if ( argc > 1 && strcmp( argv[1], "vegas" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::Type::Vegas; if ( argc > 1 && strcmp( argv[1], "miser" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::Type::MISER; { cout << "Testing with " << mg.parameters->integrator.type << " integrator" << endl; } mg.parameters->kinematics.setSqrtS( 13.e3 ); mg.parameters->kinematics.cuts.central.eta_single.in( -2.5, 2.5 ); mg.parameters->kinematics.cuts.remnants.mass_single.max() = 1000.; //mg.parameters->integrator.ncvg = 50000; CG_INFO( "main" ) << "Initial configuration time: " << tmr.elapsed()*1.e3 << " ms."; tmr.reset(); unsigned short num_tests = 0, num_tests_passed = 0; vector<string> failed_tests, passed_tests; try { for ( const auto& values_vs_generator : values_map ) { // loop over all generators const string generator = values_vs_generator.first; if ( generator == "lpair" ) mg.parameters->setProcess( new CepGen::Process::GamGamLL ); else if ( generator == "pptoll" ) { mg.parameters->setProcess( new CepGen::Process::PPtoFF ); mg.parameters->kinematics.central_system = { CepGen::PDG::Muon, CepGen::PDG::Muon }; mg.parameters->kinematics.cuts.initial.qt = { 0., 50. }; } else if ( generator == "pptoww" ) { mg.parameters->setProcess( new CepGen::Process::PPtoWW ); mg.parameters->kinematics.setSqrtS( 13.e3 ); //mg.parameters->kinematics.cuts.initial.qt = { 0., 50. }; } else { CG_ERROR( "main" ) << "Unrecognized generator mode: " << values_vs_generator.first << "."; break; } for ( const auto& values_vs_cut : values_vs_generator.second ) { // loop over the single lepton pT cut mg.parameters->kinematics.cuts.central.pt_single.min() = values_vs_cut.first; for ( const auto& values_vs_kin : values_vs_cut.second ) { // loop over all possible kinematics const string kin_mode = values_vs_kin.first; if ( kin_mode.find( "elastic" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::Mode::ElasticElastic; else if ( kin_mode.find( "singlediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::Mode::InelasticElastic; else if ( kin_mode.find( "doublediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::Mode::InelasticInelastic; else { CG_ERROR( "main" ) << "Unrecognized kinematics mode: " << values_vs_kin.first << "."; break; } if ( kin_mode.find( "_su" ) != string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SzczurekUleshchenko; else if ( kin_mode.find( "_lux" ) != string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::Schaefer; else if ( kin_mode.find( "_allm" ) != string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::ALLM97; else mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SuriYennie; //mg.parameters->dump(); CG_INFO( "main" ) << "Process: "<< values_vs_generator.first << "/" << values_vs_kin.first << "\n\t" << "Configuration time: " << tmr.elapsed()*1.e3 << " ms."; tmr.reset(); mg.clearRun(); const double xsec_ref = values_vs_kin.second.first; const double err_xsec_ref = values_vs_kin.second.second; double xsec_cepgen, err_xsec_cepgen; mg.computeXsection( xsec_cepgen, err_xsec_cepgen ); const double sigma = fabs( xsec_ref-xsec_cepgen ) / std::hypot( err_xsec_cepgen, err_xsec_ref ); CG_INFO( "main" ) << "Computed cross section:\n\t" << "Ref. = " << xsec_ref << " +/- " << err_xsec_ref << "\n\t" << "CepGen = " << xsec_cepgen << " +/- " << err_xsec_cepgen << "\n\t" << "Pull: " << sigma << "."; CG_INFO( "main" ) << "Computation time: " << tmr.elapsed()*1.e3 << " ms."; tmr.reset(); ostringstream oss; oss << values_vs_kin.first; string test_res = Form( "%-10s", values_vs_generator.first )+"\t"+ Form( "pt-gt-%.1f", values_vs_cut.first )+"\t"+ Form( "%-16s", oss.str().c_str() )+"\t" "ref="+Form( "%g", xsec_ref )+"\t" "got="+Form( "%g", xsec_cepgen )+"\t" "pull="+Form( "%+g", sigma ); if ( fabs( sigma ) < num_sigma ) { passed_tests.emplace_back( test_res ); num_tests_passed++; } else failed_tests.emplace_back( test_res ); num_tests++; cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << endl; } } } } catch ( CepGen::Exception& e ) {} if ( failed_tests.size() != 0 ) { ostringstream os_failed, os_passed; for ( const auto& fail : failed_tests ) os_failed << " (*) " << fail << endl; for ( const auto& pass : passed_tests ) os_passed << " (*) " << pass << endl; throw CG_FATAL( "main" ) << "Some tests failed!\n" << os_failed.str() << "\n" << "Passed tests:\n" << os_passed.str() << "."; } CG_INFO( "main" ) << "ALL TESTS PASSED!"; return 0; } diff --git a/test/test_singlepoint.cpp b/test/test_singlepoint.cpp index cd4864a..00495a3 100644 --- a/test/test_singlepoint.cpp +++ b/test/test_singlepoint.cpp @@ -1,31 +1,31 @@ #include "CepGen/Generator.h" #include "CepGen/Parameters.h" #include "CepGen/Processes/GamGamLL.h" using namespace std; int main() { CepGen::Generator g; CepGen::Parameters* p = g.parameters.get(); //p->setProcess( new GamGamLL ); p->setProcess( new CepGen::Process::GamGamLL ); p->kinematics.mode = CepGen::Kinematics::Mode::ElasticElastic; //p->kinematics.mode = CepGen::Kinematics::Mode::InelasticElastic; //p->kinematics.mode = CepGen::Kinematics::Mode::ElasticInelastic; p->kinematics.cuts.central[CepGen::Cuts::pt_single] = 5.; p->kinematics.cuts.central[CepGen::Cuts::eta_single] = { -2.5, 2.5 }; p->kinematics.cuts.remnants[CepGen::Cuts::mass_single] = { 1.07, 320. }; p->dump(); - CepGen::Logger::get().level = CepGen::Logger::Level::DebugInsideLoop; + CepGen::Logger::get().level = CepGen::Logger::Level::debugInsideLoop; const unsigned short ndim = g.numDimensions(); double x[12]; for ( unsigned int i = 0; i < ndim; ++i ) x[i] = 0.3; cout << g.computePoint( x ) << endl; return 0; }