diff --git a/Analysis/EventShapes.icc b/Analysis/EventShapes.icc --- a/Analysis/EventShapes.icc +++ b/Analysis/EventShapes.icc @@ -1,317 +1,316 @@ // -*- C++ -*- // // This is the implementation of the inlined member functions of // the EventShapes class. // namespace ThePEG { inline EventShapes::EventShapes() {} inline EventShapes::EventShapes(const EventShapes & x) : Interfaced(x) {} inline IBPtr EventShapes::clone() const { return new_ptr(*this); } inline IBPtr EventShapes::fullclone() const { return new_ptr(*this); } -inline void EventShapes::doupdate() throw(UpdateException) { +inline void EventShapes::doupdate() { Interfaced::doupdate(); // First update base class. bool redo = touched(); // redo if touched. // UpdateChecker::check(aDependentMember, redo); // Update referenced objects on which this depends redo is set to true // if the dependent object is touched. // for_each(ContainerOfDependencies, UpdateChecker(redo)); // Update a container of references. // for_each(MapOfDependencies, UpdateMapChecker(redo)); // Update a map of references. if ( !redo ) return; // return if nothing has been touched. Otherwise do the actual update. // touch() // Touch if anything has changed. } -inline void EventShapes::doinit() throw(InitException) { +inline void EventShapes::doinit() { Interfaced::doinit(); } inline void EventShapes::dofinish() { Interfaced::dofinish(); } inline void EventShapes::doinitrun() { Interfaced::doinitrun(); } -inline void EventShapes::rebind(const TranslationMap & trans) - throw(RebindException) { +inline void EventShapes::rebind(const TranslationMap & trans) { // dummy = trans.translate(dummy); Interfaced::rebind(trans); } inline IVector EventShapes::getReferences() { IVector ret = Interfaced::getReferences(); // ret.push_back(dummy); return ret; } inline double EventShapes::thrust() { checkThrust(); return _thrust[0]; } inline double EventShapes::thrustMajor() { checkThrust(); return _thrust[1]; } inline double EventShapes::thrustMinor() { checkThrust(); return _thrust[2]; } inline double EventShapes::oblateness() { checkThrust(); return _thrust[1]-_thrust[2]; } inline ThreeVector EventShapes::thrustAxis() { checkThrust(); return _thrustAxis[0]; } inline ThreeVector EventShapes::majorAxis() { checkThrust(); return _thrustAxis[1]; } inline ThreeVector EventShapes::minorAxis() { checkThrust(); return _thrustAxis[2]; } inline void EventShapes::reset(const tPVector &part) { _pv.resize(part.size()); for(unsigned int ix=0;ixmomentum(); _thrustDone = false; _spherDone = false; _linTenDone = false; _hemDone = false; _useCmBoost = false; } inline double EventShapes::sphericity() { checkSphericity(); return 3./2.*(_spher[1]+_spher[2]); } inline double EventShapes::aplanarity() { checkSphericity(); return 3./2.*_spher[2]; } inline double EventShapes::planarity() { checkSphericity(); return _spher[1]-_spher[2]; } inline ThreeVector EventShapes::sphericityAxis() { checkSphericity(); return _spherAxis[0]; } inline vector EventShapes::sphericityEigenValues() { checkSphericity(); return _spher; } inline vector< ThreeVector > EventShapes::sphericityEigenVectors() { checkSphericity(); return _spherAxis; } inline vector EventShapes::linTenEigenValues() { checkLinTen(); return _linTen; } inline vector< ThreeVector > EventShapes::linTenEigenVectors() { checkLinTen(); return _linTenAxis; } inline double EventShapes::CParameter() { checkLinTen(); return 3.*(_linTen[0]*_linTen[1]+_linTen[1]*_linTen[2] +_linTen[2]*_linTen[0]); } inline double EventShapes::DParameter() { checkLinTen(); return 27.*(_linTen[0]*_linTen[1]*_linTen[2]); } inline double EventShapes::Mhigh2() { checkHemispheres(); return _mPlus; } inline double EventShapes::Mlow2() { checkHemispheres(); return _mMinus; } inline double EventShapes::Mdiff2() { checkHemispheres(); return _mPlus-_mMinus; } inline double EventShapes::Bmax() { checkHemispheres(); return _bPlus; } inline double EventShapes::Bmin() { checkHemispheres(); return _bMinus; } inline double EventShapes::Bsum() { checkHemispheres(); return _bPlus+_bMinus; } inline double EventShapes::Bdiff() { checkHemispheres(); return _bPlus-_bMinus; } inline void EventShapes::checkLinTen() { if (!_linTenDone) { _linTenDone = true; diagonalizeTensors(true, _useCmBoost); } } inline void EventShapes::checkSphericity() { if (!_spherDone) { _spherDone = true; diagonalizeTensors(false, _useCmBoost); } } inline void EventShapes::checkThrust() { if (!_thrustDone) { _thrustDone = true; calculateThrust(); } } inline void EventShapes::checkHemispheres() { if (!_hemDone) { _hemDone = true; calcHemisphereMasses(); } } inline void EventShapes::calcHemisphereMasses() { Lorentz5Momentum pos, neg; Energy pden(0.0*GeV), epos(0.0*GeV), eneg(0.0*GeV); for(unsigned int ix=0;ix<_pv.size();++ix) { if( _pv[ix].vect()*thrustAxis() > 0.0*GeV) { pos += _pv[ix]; epos += _pv[ix].perp(thrustAxis()); } else { neg += _pv[ix]; eneg += _pv[ix].perp(thrustAxis()); } pden += _pv[ix].vect().mag(); } // denominator and masses Energy2 den(sqr(pos.e()+neg.e())); _mPlus = pos.m2()/den; _mMinus = neg.m2()/den; if (_mPlus < _mMinus) swap(_mPlus, _mMinus); // jet broadening _bPlus = 0.5*epos/pden; _bMinus = 0.5*eneg/pden; if (_bPlus < _bMinus) swap(_bPlus, _bMinus); } inline double EventShapes::getX(const Lorentz5Momentum & p, const Energy & Ebeam) { return ( Ebeam > 0.0*GeV? double(p.vect().mag()/Ebeam): -1.0 ); } inline double EventShapes::getXi(const Lorentz5Momentum & p, const Energy & Ebeam) { return( ( Ebeam > 0.0*GeV && p.vect().mag() > 0.0*GeV)? log(Ebeam/p.vect().mag()) : -1.0 ); } inline Energy EventShapes::getPt(const Lorentz5Momentum & p) { return p.perp(); } inline double EventShapes::getRapidity(const Lorentz5Momentum & p) { return (p.t() > p.z() ? p.rapidity() : 1e99); } inline Energy EventShapes::ptInT(const Lorentz5Momentum & p) { checkThrust(); return p.vect()*_thrustAxis[1]; } inline Energy EventShapes::ptOutT(const Lorentz5Momentum & p) { checkThrust(); return p.vect()*_thrustAxis[2]; } inline double EventShapes::yT(const Lorentz5Momentum & p) { checkThrust(); return (p.t() > p.vect()*_thrustAxis[0] ? p.rapidity(_thrustAxis[0]) : 1e99); } inline Energy EventShapes::ptInS(const Lorentz5Momentum & p) { checkSphericity(); return p.vect()*_spherAxis[1]; } inline Energy EventShapes::ptOutS(const Lorentz5Momentum & p) { checkSphericity(); return p.vect()*_spherAxis[2]; } inline double EventShapes::yS(const Lorentz5Momentum & p) { checkSphericity(); return (p.t() > p.vect()*_spherAxis[0] ? p.rapidity(_spherAxis[0]) : 1e99); } inline void EventShapes::normalizeEEC(vector & hi, long evts) { for (unsigned int bin = 0; bin < hi.size(); bin++) bin /= (hi.size()*evts); } inline double EventShapes::AEEC(vector & hi, double& coschi) { if (coschi > 0. && coschi <= 1.) { int i = (int) floor((-coschi+1.)/2.*hi.size()); int j = (int) floor((coschi+1.)/2.*hi.size()); return hi[i]-hi[j]; } else { return 1e99; } } } diff --git a/Analysis/LEPEventShapes.cc b/Analysis/LEPEventShapes.cc --- a/Analysis/LEPEventShapes.cc +++ b/Analysis/LEPEventShapes.cc @@ -1,538 +1,538 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the LEPEventShapes class. // #include "LEPEventShapes.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Repository/CurrentGenerator.h" #ifdef ThePEG_TEMPLATES_IN_CC_FILE // #include "LEPEventShapes.tcc" #endif #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #ifndef LWH_AIAnalysisFactory_H #ifndef LWH #define LWH ThePEGLWH #endif #include "ThePEG/Analysis/LWH/AnalysisFactory.h" #endif using namespace ThePEG; using namespace ThePEG; void LEPEventShapes::analyze(tEventPtr event, long ieve, int loop, int state) { AnalysisHandler::analyze(event, ieve, loop, state); if ( loop > 0 || state != 0 || !event ) return; // get the final-state particles tPVector hadrons=event->getFinalState(); // event shapes } LorentzRotation LEPEventShapes::transform(tEventPtr event) const { return LorentzRotation(); // Return the Rotation to the frame in which you want to perform the analysis. } void LEPEventShapes::analyze(const tPVector & particles) { double eventweight = generator()->currentEvent()->weight(); _shapes->reset(particles); _omthr ->fill( 1.-_shapes->thrust() ,eventweight); _maj ->fill( _shapes->thrustMajor() ,eventweight); _min ->fill( _shapes->thrustMinor() ,eventweight); _obl ->fill( _shapes->oblateness() ,eventweight); _c ->fill( _shapes->CParameter() ,eventweight); _d ->fill( _shapes->DParameter() ,eventweight); _sph ->fill( _shapes->sphericity() ,eventweight); _apl ->fill( _shapes->aplanarity() ,eventweight); _pla ->fill( _shapes->planarity() ,eventweight); _mhi ->fill( _shapes->Mhigh2() ,eventweight); _mlo ->fill( _shapes->Mlow2() ,eventweight); _mdiff ->fill( _shapes->Mdiff2() ,eventweight); _bmax ->fill( _shapes->Bmax() ,eventweight); _bmin ->fill( _shapes->Bmin() ,eventweight); _bsum ->fill( _shapes->Bsum() ,eventweight); _bdiff ->fill( _shapes->Bdiff() ,eventweight); } void LEPEventShapes::dofinish() { AnalysisHandler::dofinish(); if ( !checkHistogramFactory() ) return; unitNormalize(_omthr); unitNormalize(_maj); unitNormalize(_min); unitNormalize(_obl); unitNormalize(_sph); unitNormalize(_apl); unitNormalize(_pla); unitNormalize(_c); unitNormalize(_d); unitNormalize(_mhi); unitNormalize(_mlo); unitNormalize(_mdiff); unitNormalize(_bmax); unitNormalize(_bmin); unitNormalize(_bsum); unitNormalize(_bdiff); } void LEPEventShapes::dataSet(const std::string & path, const std::string & title, const std::vector & bins, const std::vector & y, const std::vector & yerr){ std::vector x; for(std::vector::const_iterator it = bins.begin(); it + 1 != bins.end(); it++){ x.push_back((*it + *(it+1)) / 2.0); } std::vector xerr(x.size(), 0.0); histogramFactory().dataSetFactory().createXY(path, title, x, y, xerr, yerr); } void LEPEventShapes::doinitrun() { AnalysisHandler::doinitrun(); histogramFactory().registerClient(this); histogramFactory().mkdirs("/LEPEventShapes"); histogramFactory().cd("/LEPEventShapes"); vector bins,data,error; // 1-T double vals1[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090, 0.100, 0.120, 0.140, 0.160, 0.180, 0.200, 0.250, 0.300, 0.350, 0.400, 0.500}; double data1[]= { 1.030,10.951,17.645,14.192,10.009, 7.572, 5.760, 4.619, 3.792, 3.176, 2.456, 1.825, 1.401, 1.074, 0.8262, 0.5525,0.3030,0.1312,0.0238,0.0007}; double error1stat[]={0.019 ,0.051 ,0.066 ,0.061 ,0.050 , 0.044 ,0.038 ,0.034 ,0.031 ,0.028 , 0.018 ,0.015 ,0.013 ,0.011 ,0.0100 , 0.0051 ,0.0038 ,0.0025 ,0.0012 ,0.0002 }; double error1syst[]={0.076 , 0.527 , 0.547 , 0.292 , 0.152 , 0.101 , 0.076 , 0.062 , 0.051 , 0.042 , 0.032 , 0.022 , 0.016 , 0.011 , 0.0083, 0.0065 , 0.0058, 0.0044, 0.0014, 0.0001}; double error1[20]; for(unsigned int ix=0;ix<20;++ix){error1[ix]=sqrt(sqr(error1stat[ix])+ sqr(error1syst[ix]));} bins = vector(vals1 ,vals1 +21); data = vector(data1 ,data1 +20); error = vector(error1,error1+20); _omthr = histogramFactory().createHistogram1D ("omthr", "One Minus Thrust", bins); dataSet("omthr_data", "DELPHI data", bins, data, error); // major double vals2[] = {0.000, 0.020, 0.040, 0.050, 0.060, 0.070, 0.080, 0.100, 0.120, 0.140, 0.160, 0.200, 0.240, 0.280, 0.320, 0.360, 0.400, 0.440, 0.480, 0.520, 0.560, 0.600, 0.640}; double data2[]={0.00040 ,0.0590 ,0.642 ,2.178 ,4.303 , 5.849 ,6.889 ,6.342 ,4.890 ,3.900 , 2.960 ,2.124 ,1.5562 ,1.1807 ,0.8693, 0.6493 ,0.4820 ,0.3493 ,0.2497 ,0.1489, 0.0714 ,0.0203}; double error2stat[]={0.00090 ,0.0030 ,0.013 ,0.024 ,0.034 , 0.039 ,0.030 ,0.028 ,0.024 ,0.021 , 0.013 ,0.011 ,0.0095 ,0.0083 ,0.0071 , 0.0061 ,0.0052 ,0.0044 ,0.0037 ,0.0028 , 0.0019 ,0.0010}; double error2syst[]={0.00005 ,0.0058 ,0.028 ,0.086 ,0.155 , 0.192 ,0.194 ,0.143 ,0.085 ,0.050 , 0.030 ,0.021 ,0.0156 ,0.0118 ,0.0087 , 0.0065 ,0.0048 ,0.0055 ,0.0065 ,0.0058 , 0.0038 ,0.0014}; double error2[22]; for(unsigned int ix=0;ix<22;++ix){error2[ix]=sqrt(sqr(error2stat[ix])+ sqr(error2syst[ix]));} bins = vector(vals2 ,vals2 +23); data = vector(data2 ,data2 +22); error = vector(error2,error2+22); _maj = histogramFactory().createHistogram1D ("maj", "Thrust Major", bins); dataSet("maj_data", "DELPHI data", bins, data, error); // minor double vals3[] = {0.000, 0.020, 0.040, 0.050, 0.060, 0.070, 0.080, 0.100, 0.120, 0.140, 0.160, 0.200, 0.240, 0.280, 0.320, 0.400}; double data3[]={ 0.0156 , 1.236 , 5.706 , 9.714 ,12.015 , 12.437 ,10.404 , 6.918 , 4.250 , 2.517 , 1.2561 , 0.4895 , 0.2112 , 0.0879 , 0.0250 }; double error3stat[]={0.0017 ,0.013 ,0.037 ,0.048 ,0.054 , 0.055 ,0.036 ,0.029 ,0.023 ,0.017 , 0.0086 ,0.0054 ,0.0036 ,0.0023 ,0.0009}; double error3syst[]={0.0036,0.066 ,0.073 ,0.125 ,0.155 , 0.161 ,0.136 ,0.092 ,0.058 ,0.035 , 0.0187,0.0080,0.0039,0.0018,0.0006}; double error3[15]; for(unsigned int ix=0;ix<15;++ix){error3[ix]=sqrt(sqr(error3stat[ix])+ sqr(error3syst[ix]));} bins = vector(vals3 ,vals3 +16); data = vector(data3 ,data3 +15); error = vector(error3,error3+15); _min = histogramFactory().createHistogram1D ("min", "Thrust Minor", bins); dataSet("min_data", "DELPHI data", bins, data, error); // oblateness double vals4[] = {0.000, 0.020, 0.040, 0.060, 0.080, 0.100, 0.120, 0.140, 0.160, 0.200, 0.240, 0.280, 0.320, 0.360, 0.400, 0.440, 0.520}; double data4[]={ 9.357 ,11.508 , 7.215 , 4.736 , 3.477 , 2.696 , 2.106 , 1.690 , 1.2648 , 0.8403 , 0.5674 , 0.3842 , 0.2573 , 0.1594 , 0.0836 , 0.0221 }; double error4stat[]={0.036 ,0.038 ,0.029 ,0.023 ,0.020 , 0.018 ,0.016 ,0.014 ,0.0085 ,0.0069 , 0.0056 ,0.0046 ,0.0037 ,0.0029 ,0.0020 , 0.0007}; double error4syst[]={0.178 ,0.140 ,0.072 ,0.047 ,0.035 , 0.027 ,0.021 ,0.017 ,0.0126,0.0087, 0.0065,0.0050,0.0043,0.0037,0.0030, 0.0015}; double error4[16]; for(unsigned int ix=0;ix<16;++ix){error4[ix]=sqrt(sqr(error4stat[ix])+ sqr(error4syst[ix]));} bins = vector(vals4 ,vals4 +17); data = vector(data4 ,data4 +16); error = vector(error4,error4+16); _obl = histogramFactory().createHistogram1D ("obl", "Oblateness", bins); dataSet("obl_data", "DELPHI data", bins, data, error); // sphericity double vals5[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.080, 0.100, 0.120, 0.160, 0.200, 0.250, 0.300, 0.350, 0.400, 0.500, 0.600, 0.700, 0.850}; double data5[]={16.198 ,20.008 ,12.896 , 8.237 , 5.885 , 4.458 , 3.272 , 2.290 , 1.699 , 1.2018 , 0.7988 , 0.5610 , 0.3926 , 0.2810 , 0.2099 , 0.1441 , 0.0842 , 0.04160 , 0.00758 }; double error5stat[]={0.067 ,0.072 ,0.056 ,0.043 ,0.037 , 0.032 ,0.019 ,0.016 ,0.014 ,0.0082 , 0.0067 ,0.0050 ,0.0042 ,0.0035 ,0.0030 , 0.0018 ,0.0013 ,0.00092 ,0.00032}; double error5syst[]={0.208 ,0.246 ,0.153 ,0.094 ,0.065 , 0.048 ,0.034 ,0.023 ,0.017 ,0.0120 , 0.0080 ,0.0063 ,0.0051 ,0.0043 ,0.0037 , 0.0032 ,0.0023 ,0.00129,0.00024}; double error5[19]; for(unsigned int ix=0;ix<19;++ix){error5[ix]=sqrt(sqr(error5stat[ix])+ sqr(error5syst[ix]));} bins = vector(vals5 ,vals5 +20); data = vector(data5 ,data5 +19); error = vector(error5,error5+19); _sph = histogramFactory().createHistogram1D ("sph", "Sphericity", bins); dataSet("sph_data", "DELPHI data", bins, data, error); // aplanarity double vals6[] = {0.000, 0.005, 0.010, 0.015, 0.020, 0.030, 0.040, 0.060, 0.080, 0.100, 0.120, 0.140, 0.160, 0.200, 0.250, 0.300}; double data6[]={75.10 ,55.31 ,26.03 ,13.927 , 6.768 , 3.014 , 1.281 , 0.5181 , 0.2619 , 0.1461 , 0.0758 , 0.0467 , 0.0234 , 0.00884 , 0.00310 }; double error6stat[]={0.19 ,0.17 ,0.11 ,0.079 ,0.038 , 0.025 ,0.012 ,0.0075 ,0.0054 ,0.0041 , 0.0029 ,0.0023 ,0.0011 ,0.00061 ,0.00040 }; double error6syst[]={0.75 ,0.55 ,0.28 ,0.176 ,0.098 , 0.056 ,0.035 ,0.0188 ,0.0118 ,0.0079 , 0.0043 ,0.0027 ,0.0014 ,0.00052,0.00018}; double error6[19]; - for(unsigned int ix=0;ix<19;++ix){error6[ix]=sqrt(sqr(error6stat[ix])+ + for(unsigned int ix=0;ix<15;++ix){error6[ix]=sqrt(sqr(error6stat[ix])+ sqr(error6syst[ix]));} bins = vector(vals6 ,vals6 +16); data = vector(data6 ,data6 +15); error = vector(error6,error6+15); _apl = histogramFactory().createHistogram1D ("apl", "Aplanarity", bins); dataSet("apl_data", "DELPHI data", bins, data, error); // planarity double vals7[] = {0.000, 0.005, 0.010, 0.015, 0.020, 0.025, 0.030, 0.035, 0.040, 0.050, 0.060, 0.080, 0.100, 0.120, 0.160, 0.200, 0.250, 0.300, 0.350, 0.400, 0.500}; double data7[]={68.69 ,31.66 ,17.091 ,11.370 , 8.417 , 6.578 , 5.479 , 4.493 , 3.610 , 2.749 , 1.987 , 1.362 , 1.008 , 0.6676 , 0.4248 , 0.2692 , 0.1742 , 0.1042 , 0.0566 , 0.0145 }; double error7stat[]={0.19 ,0.12 ,0.088 ,0.072 ,0.062 , 0.055 ,0.050 ,0.045 ,0.029 ,0.025 , 0.015 ,0.012 ,0.011 ,0.0061 ,0.0048 , 0.0034 ,0.0028 ,0.0021 ,0.0015 ,0.0006}; double error7syst[]={0.74 ,0.35 ,0.188 ,0.127 ,0.095 , 0.075 ,0.063 ,0.052 ,0.042 ,0.033 , 0.024 ,0.017 ,0.013 ,0.0093,0.0063, 0.0042 ,0.0029,0.0019,0.0011,0.0003}; double error7[20]; for(unsigned int ix=0;ix<20;++ix){error7[ix]=sqrt(sqr(error7stat[ix])+ sqr(error7syst[ix]));} bins = vector(vals7 ,vals7 +21); data = vector(data7 ,data7 +20); error = vector(error7,error7+20); _pla = histogramFactory().createHistogram1D ("pla", "Planarity", bins); dataSet("pla_data", "DELPHI data", bins, data, error); // C double vals8[] = {0.000, 0.040, 0.080, 0.120, 0.160, 0.200, 0.240, 0.280, 0.320, 0.360, 0.400, 0.440, 0.480, 0.520, 0.560, 0.600, 0.640, 0.680, 0.720, 0.760, 0.800, 0.840, 0.880, 0.920}; double data8[]={0.0881 ,1.5383 ,3.909 ,3.833 ,2.835 , 2.164 ,1.716 ,1.3860 ,1.1623 ,0.9720 , 0.8349 ,0.7161 ,0.6205 ,0.5441 ,0.4844 , 0.4209 ,0.3699 ,0.3286 ,0.2813 ,0.2178 , 0.1287 ,0.0542 ,0.0212 }; double error8stat[]={0.0030 ,0.0100 ,0.016 ,0.016 ,0.013 , 0.012 ,0.010 ,0.0092 ,0.0084 ,0.0077 , 0.0072 ,0.0066 ,0.0061 ,0.0057 ,0.0054 , 0.0050 ,0.0046 ,0.0044 ,0.0040 ,0.0033 , 0.0026 ,0.0016 ,0.0009}; double error8syst[]={0.0067,0.0831,0.142 ,0.088 ,0.040 , 0.022 ,0.017 ,0.0139,0.0116,0.0097, 0.0083,0.0072,0.0062,0.0054,0.0050, 0.0063,0.0079,0.0099,0.0129,0.0151, 0.0130,0.0076,0.0040}; double error8[23]; for(unsigned int ix=0;ix<23;++ix){error8[ix]=sqrt(sqr(error8stat[ix])+ sqr(error8syst[ix]));} bins = vector(vals8 ,vals8 +24); data = vector(data8 ,data8 +23); error = vector(error8,error8+23); _c = histogramFactory().createHistogram1D ("c", "C parameter", bins); dataSet("c_data", "DELPHI data", bins, data, error); // D double vals9[] = {0.000, 0.008, 0.016, 0.030, 0.044, 0.066, 0.088, 0.112, 0.136, 0.162, 0.188, 0.218, 0.248, 0.284, 0.320, 0.360, 0.400, 0.450, 0.500, 0.560, 0.620, 0.710, 0.800}; double data9[]={22.228 ,22.766 ,12.107 , 6.879 , 4.284 , 2.727 , 1.909 , 1.415 , 1.051 , 0.7977 , 0.6155 , 0.4566 , 0.3341 , 0.2452 , 0.1774 , 0.1234 , 0.0902 , 0.0603 , 0.0368 , 0.0222 , 0.0128 , 0.0052 }; double error9stat[]={0.082 ,0.085 ,0.047 ,0.035 ,0.022 , 0.018 ,0.014 ,0.012 ,0.010 ,0.0089 , 0.0073 ,0.0063 ,0.0049 ,0.0042 ,0.0033 , 0.0028 ,0.0021 ,0.0017 ,0.0012 ,0.0009 , 0.0006 ,0.0004}; double error9syst[]={0.868 ,0.440 ,0.150 ,0.079 ,0.053 , 0.036 ,0.028 ,0.022 ,0.018 ,0.0145, 0.0117 ,0.0089,0.0065,0.0049,0.0037 , 0.0028,0.0023 ,0.0018,0.0013,0.0009, 0.0006,0.0003}; double error9[22]; for(unsigned int ix=0;ix<22;++ix){error9[ix]=sqrt(sqr(error9stat[ix])+ sqr(error9syst[ix]));} bins = vector(vals9 ,vals9 +23); data = vector(data9 ,data9 +22); error = vector(error9,error9+22); _d = histogramFactory().createHistogram1D ("d", "D parameter", bins); dataSet("d_data", "DELPHI data", bins, data, error); // M high double vals10[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.080, 0.100, 0.120, 0.140, 0.160, 0.200, 0.250, 0.300, 0.350, 0.400}; double data10[]={ 1.994 ,18.580 ,20.678 ,13.377 , 8.965 , 6.558 , 4.515 , 2.914 , 1.991 , 1.406 , 1.010 , 0.6319 , 0.3085 , 0.1115 , 0.0184 , 0.0008 }; double error10stat[]={0.027 ,0.065 ,0.076 ,0.060 ,0.049 , 0.041 ,0.024 ,0.019 ,0.016 ,0.013 , 0.011 ,0.0063 ,0.0039 ,0.0022 ,0.0008 , 0.0002 }; double error10syst[]={0.166 ,0.709 ,0.729 ,0.412 ,0.239 , 0.151 ,0.082 ,0.037 ,0.020 ,0.014 , 0.010 ,0.0063,0.0051,0.0039,0.0012, 0.0001}; double error10[16]; for(unsigned int ix=0;ix<16;++ix){error10[ix]=sqrt(sqr(error10stat[ix])+ sqr(error10syst[ix]));} bins = vector(vals10 ,vals10 +17); data = vector(data10 ,data10 +16); error = vector(error10,error10+16); _mhi = histogramFactory().createHistogram1D ("mhi", "High hemisphere mass", bins); dataSet("mhi_data", "DELPHI data", bins, data, error); // M low double vals11[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.080, 0.100, 0.120}; double data11[]={23.414 ,39.12 ,18.080 , 7.704 , 3.922 , 2.128 , 1.013 , 0.3748 , 0.1412 }; double error11stat[]={0.074 ,0.11 ,0.081 ,0.052 ,0.036 , 0.026 ,0.013 ,0.0079 ,0.0050}; double error11syst[]={1.595 ,2.65 ,1.215 ,0.514 ,0.260 , 0.140 ,0.066 ,0.0241,0.0089}; double error11[9]; for(unsigned int ix=0;ix<9;++ix){error11[ix]=sqrt(sqr(error11stat[ix])+ sqr(error11syst[ix]));} bins = vector(vals11 ,vals11 +10); data = vector(data11 ,data11 + 9); error = vector(error11,error11+ 9); _mlo = histogramFactory().createHistogram1D ("mlo", "Low hemisphere mass", bins); dataSet("mlo_data", "DELPHI data", bins, data, error); // M diff double vals12[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.060, 0.080, 0.120, 0.160, 0.200, 0.250, 0.300, 0.350, 0.400}; double data12[]={35.393 ,20.745 ,11.426 , 7.170 , 4.344 , 2.605 , 1.4238 , 0.7061 , 0.3831 , 0.1836 , 0.0579 , 0.0075 , 0.0003}; double error12stat[]={0.092 ,0.071 ,0.052 ,0.041 ,0.023 , 0.017 ,0.0092 ,0.0064 ,0.0046 ,0.0028 , 0.0015 ,0.0006 ,0.0002}; double error12syst[]={0.354 ,0.207 ,0.114 ,0.072 ,0.043 , 0.026 ,0.0142,0.0071,0.0044,0.0032, 0.0018,0.0006,0.0001}; double error12[13]; for(unsigned int ix=0;ix<13;++ix){error12[ix]=sqrt(sqr(error12stat[ix])+ sqr(error12syst[ix]));} bins = vector(vals12 ,vals12 +14); data = vector(data12 ,data12 +13); error = vector(error12,error12+13); //_mdiff= new_ptr(Histogram(bins,data,error)); _mdiff = histogramFactory().createHistogram1D ("mdiff", "Difference in hemisphere masses", bins); dataSet("mdiff_data", "DELPHI data", bins, data, error); // Bmax double vals13[] = {0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080, 0.100, 0.120, 0.140, 0.170, 0.200, 0.240, 0.280, 0.320}; double data13[]={0.6707 , 7.538 ,14.690 ,13.942 ,11.298 , 9.065 , 7.387 , 5.445 , 3.796 , 2.670 , 1.756 , 1.0580 , 0.5288 , 0.1460 , 0.0029 }; double error13stat[]={0.0096 ,0.038 ,0.058 ,0.057 ,0.053 , 0.048 ,0.043 ,0.026 ,0.022 ,0.018 , 0.012 ,0.0092 ,0.0056 ,0.0028 ,0.0004}; double error13syst[]={0.1077,0.809 ,0.745 ,0.592 ,0.379 , 0.266 ,0.222 ,0.176 ,0.127 ,0.087 , 0.051 ,0.0218,0.0053,0.0071,0.0003}; double error13[15]; for(unsigned int ix=0;ix<15;++ix){error13[ix]=sqrt(sqr(error13stat[ix])+ sqr(error13syst[ix]));} bins = vector(vals13 ,vals13 +16); data = vector(data13 ,data13 +15); error = vector(error13,error13+15); _bmax = histogramFactory().createHistogram1D ("bmax", "Wide jet broadening measure", bins); dataSet("bmax_data", "DELPHI data", bins, data, error); // Bmin double vals14[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.080, 0.100, 0.120, 0.150, 0.180}; double data14[]={0.645 ,11.169 ,28.908 ,25.972 ,14.119 , 7.500 , 3.405 , 1.320 , 0.5448 , 0.1916 , 0.0366}; double error14stat[]={0.010 ,0.045 ,0.082 ,0.083 ,0.061 , 0.044 ,0.021 ,0.013 ,0.0082 ,0.0040 , 0.0017}; double error14syst[]={0.096 ,1.006 ,1.823 ,1.478 ,0.860 , 0.494 ,0.233 ,0.089 ,0.0328,0.0104, 0.0034}; double error14[11]; for(unsigned int ix=0;ix<11;++ix){error14[ix]=sqrt(sqr(error14stat[ix])+ sqr(error14syst[ix]));} bins = vector(vals14 ,vals14 +12); data = vector(data14 ,data14 +11); error = vector(error14,error14+11); _bmin = histogramFactory().createHistogram1D ("bmin", "Narrow jet broadening measure", bins); dataSet("bmin_data", "DELPHI data", bins, data, error); // Bsum double vals15[] = {0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090, 0.100, 0.110, 0.130, 0.150, 0.170, 0.190, 0.210, 0.240, 0.270, 0.300, 0.330, 0.360}; double data15[]={0.2030 ,1.628 ,4.999 ,8.190 ,9.887 , 9.883 ,9.007 ,7.746 ,6.714 ,5.393 , 3.998 ,2.980 ,2.294 ,1.747 ,1.242 , 0.8125 ,0.4974 ,0.2285 ,0.0732 }; double error15stat[]={0.0055 ,0.015 ,0.031 ,0.041 ,0.047 , 0.049 ,0.047 ,0.044 ,0.041 ,0.026 , 0.023 ,0.019 ,0.017 ,0.015 ,0.010 , 0.0080 ,0.0062 ,0.0041 ,0.0024}; double error15syst[]={0.0383,0.183 ,0.463 ,0.644 ,0.661 , 0.564 ,0.443 ,0.332 ,0.255 ,0.180 , 0.125 ,0.098 ,0.085 ,0.075 ,0.063 , 0.0469,0.0296,0.0119,0.0007}; double error15[19]; for(unsigned int ix=0;ix<19;++ix){error15[ix]=sqrt(sqr(error15stat[ix])+ sqr(error15syst[ix]));} bins = vector(vals15 ,vals15 +20); data = vector(data15 ,data15 +19); error = vector(error15,error15+19); _bsum = histogramFactory().createHistogram1D ("bsum", "Sum of jet broadening measures", bins); dataSet("bsum_data", "DELPHI data", bins, data, error); // Bdiff double vals16[] = {0.000, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090, 0.100, 0.120, 0.140, 0.160, 0.180, 0.200, 0.240, 0.280}; double data16[]={26.630 ,18.684 ,12.343 , 8.819 , 6.688 , 5.111 , 4.071 , 3.271 , 2.681 , 2.233 , 1.647 , 1.111 , 0.7618 , 0.5138 , 0.3167 , 0.1265 , 0.0117}; double error16stat[]={0.081 ,0.066 ,0.054 ,0.046 ,0.040 , 0.035 ,0.031 ,0.028 ,0.025 ,0.023 , 0.014 ,0.011 ,0.0095 ,0.0078 ,0.0062 , 0.0026 ,0.0008}; double error16syst[]={0.459 ,0.292 ,0.186 ,0.134 ,0.106 , 0.084 ,0.068 ,0.054 ,0.043 ,0.035 , 0.026 ,0.019 ,0.0144,0.0119,0.0098, 0.0056,0.0008}; double error16[17]; for(unsigned int ix=0;ix<17;++ix){error16[ix]=sqrt(sqr(error16stat[ix])+ sqr(error16syst[ix]));} bins = vector(vals16 ,vals16 +18); data = vector(data16 ,data16 +17); error = vector(error16,error16+17); _bdiff = histogramFactory().createHistogram1D ("bdiff", "Difference of jet broadening measures", bins); dataSet("bdiff_data", "DELPHI data", bins, data, error); } void LEPEventShapes::analyze(tPPtr) {} void LEPEventShapes::persistentOutput(PersistentOStream & os) const { os << _shapes; } void LEPEventShapes::persistentInput(PersistentIStream & is, int) { is >> _shapes; } ClassDescription LEPEventShapes::initLEPEventShapes; // Definition of the static class description member. void LEPEventShapes::Init() { static ClassDocumentation documentation ("The LEPEventShapes class compares event shapes at the Z mass" "with experimental results"); static Reference interfaceEventShapes ("EventShapes", "Pointer to the object which calculates the event shapes", &LEPEventShapes::_shapes, false, false, true, false, false); } diff --git a/HerwigTest/DYMass.h b/HerwigTest/DYMass.h --- a/HerwigTest/DYMass.h +++ b/HerwigTest/DYMass.h @@ -1,238 +1,237 @@ // -*- C++ -*- #ifndef THEPEG_DYMass_H #define THEPEG_DYMass_H // // This is the declaration of the DYMass class. // #include "ThePEG/Handlers/AnalysisHandler.h" #include "DYMass.fh" namespace ThePEG { /** * Here is the documentation of the DYMass class. * * @see \ref DYMassInterfaces "The interfaces" * defined for DYMass. */ class DYMass: public AnalysisHandler { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ inline DYMass(); /** * The copy constructor. */ inline DYMass(const DYMass &); /** * The destructor. */ virtual ~DYMass(); //@} public: /** @name Virtual functions required by the AnalysisHandler class. */ //@{ /** * Analyze a given Event. Note that a fully generated event * may be presented several times, if it has been manipulated in * between. The default version of this function will call transform * to make a lorentz transformation of the whole event, then extract * all final state particles and call analyze(tPVector) of this * analysis object and those of all associated analysis objects. The * default version will not, however, do anything on events which * have not been fully generated, or have been manipulated in any * way. * @param event pointer to the Event to be analyzed. * @param ieve the event number. * @param loop the number of times this event has been presented. * If negative the event is now fully generated. * @param state a number different from zero if the event has been * manipulated in some way since it was last presented. */ virtual void analyze(tEventPtr event, long ieve, int loop, int state); /** * Transform the event to the desired Lorentz frame and return the * corresponding LorentzRotation. * @param event a pointer to the Event to be transformed. * @return the LorentzRotation used in the transformation. */ virtual LorentzRotation transform(tEventPtr event) const; /** * Analyze the given vector of particles. The default version calls * analyze(tPPtr) for each of the particles. * @param particles the vector of pointers to particles to be analyzed */ virtual void analyze(const tPVector & particles); /** * Analyze the given particle. * @param particle pointer to the particle to be analyzed. */ virtual void analyze(tPPtr particle); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ inline virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ inline virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Check sanity of the object during the setup phase. */ - inline virtual void doupdate() throw(UpdateException); + inline virtual void doupdate(); /** * Initialize this object after the setup phase before saving and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ - inline virtual void doinit() throw(InitException); + inline virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); /** * Rebind pointer to other Interfaced objects. Called in the setup phase * after all objects used in an EventGenerator has been cloned so that * the pointers will refer to the cloned objects afterwards. * @param trans a TranslationMap relating the original objects to * their respective clones. * @throws RebindException if no cloned object was found for a given * pointer. */ - inline virtual void rebind(const TranslationMap & trans) - throw(RebindException); + inline virtual void rebind(const TranslationMap & trans); /** * Return a vector of all pointers to Interfaced objects used in this * object. * @return a vector of pointers. */ inline virtual IVector getReferences(); //@} private: tH1DPtr hist; tH1DPtr vhist; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initDYMass; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DYMass & operator=(const DYMass &); }; } // CLASSDOC OFF #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DYMass. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DYMass. */ typedef AnalysisHandler NthBase; }; /** This template specialization informs ThePEG about the name of * the DYMass class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "ThePEG::DYMass"; } /** Return the name of the shared library be loaded to get * access to the DYMass class and every other class it uses * (except the base class). */ static string library() { return "DYMass.so"; } }; /** @endcond */ } #include "DYMass.icc" #ifndef ThePEG_TEMPLATES_IN_CC_FILE // #include "DYMass.tcc" #endif #endif /* THEPEG_DYMass_H */ diff --git a/HerwigTest/DYMass.icc b/HerwigTest/DYMass.icc --- a/HerwigTest/DYMass.icc +++ b/HerwigTest/DYMass.icc @@ -1,57 +1,56 @@ // -*- C++ -*- // // This is the implementation of the inlined member functions of // the DYMass class. // namespace ThePEG { inline DYMass::DYMass() : hist(0), vhist(0) {} inline DYMass::DYMass(const DYMass & x) : AnalysisHandler(x), hist(0), vhist(0) {} inline IBPtr DYMass::clone() const { return new_ptr(*this); } inline IBPtr DYMass::fullclone() const { return new_ptr(*this); } -inline void DYMass::doupdate() throw(UpdateException) { +inline void DYMass::doupdate() { AnalysisHandler::doupdate(); // First update base class. bool redo = touched(); // redo if touched. // UpdateChecker::check(aDependentMember, redo); // Update referenced objects on which this depends redo is set to true // if the dependent object is touched. // for_each(ContainerOfDependencies, UpdateChecker(redo)); // Update a container of references. // for_each(MapOfDependencies, UpdateMapChecker(redo)); // Update a map of references. if ( !redo ) return; // return if nothing has been touched. Otherwise do the actual update. // touch() // Touch if anything has changed. } -inline void DYMass::doinit() throw(InitException) { +inline void DYMass::doinit() { AnalysisHandler::doinit(); } -inline void DYMass::rebind(const TranslationMap & trans) - throw(RebindException) { +inline void DYMass::rebind(const TranslationMap & trans) { // dummy = trans.translate(dummy); AnalysisHandler::rebind(trans); } inline IVector DYMass::getReferences() { IVector ret = AnalysisHandler::getReferences(); // ret.push_back(dummy); return ret; } } diff --git a/HerwigTest/MEqq2gZ2ll.cc b/HerwigTest/MEqq2gZ2ll.cc --- a/HerwigTest/MEqq2gZ2ll.cc +++ b/HerwigTest/MEqq2gZ2ll.cc @@ -1,174 +1,174 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEqq2gZ2ll class. // #include "MEqq2gZ2ll.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace ThePEG; using namespace Herwig; MEqq2gZ2ll::MEqq2gZ2ll() : coefs(20), mZ2(0.0*GeV2), GZ2(0.0*GeV2), lastCont(0.0), lastBW(0.0) {} MEqq2gZ2ll::MEqq2gZ2ll(const MEqq2gZ2ll & x) : ME2to2QCD(x), coefs(x.coefs), mZ2(x.mZ2), GZ2(x.GZ2), lastCont(x.lastCont), lastBW(x.lastBW) {} MEqq2gZ2ll::~MEqq2gZ2ll() {} unsigned int MEqq2gZ2ll::orderInAlphaS() const { return 0; } unsigned int MEqq2gZ2ll::orderInAlphaEW() const { return 2; } void MEqq2gZ2ll::getDiagrams() const { tcPDPtr gamma = getParticleData(ParticleID::gamma); tcPDPtr Z0 = getParticleData(ParticleID::Z0); for(int i = 1; i <= maxFlavour(); ++i) { tcPDPtr q = getParticleData(i); tcPDPtr qb = q->CC(); for(long lid = ParticleID::eminus; lid <= ParticleID::tauminus; lid+=2) { tcPDPtr l = getParticleData(lid); tcPDPtr lb = l->CC(); add(new_ptr((Tree2toNDiagram(2), q, qb, 1, gamma, 3, l, 3, lb, -1))); add(new_ptr((Tree2toNDiagram(2), q, qb, 1, Z0, 3, l, 3, lb, -2))); } } } Energy2 MEqq2gZ2ll::scale() const { return sHat(); } double MEqq2gZ2ll::me2() const { Energy2 m2 = meMomenta()[2].mass2(); // Energy2 p1p3 = 0.5*(m2 - tHat()); // Energy2 p1p2 = 0.5*sHat(); Energy2 p1p3 = meMomenta()[0].dot(meMomenta()[2]); Energy2 p1p2 = meMomenta()[0].dot(meMomenta()[1]); Energy4 pt2 = sqr(p1p3); Energy4 pts = p1p3*p1p2; Energy4 ps2 = sqr(p1p2); Energy4 psm = p1p2*m2; int up = abs(mePartonData()[0]->id() + 1)%2; lastCont = (coefs[0 + up]*(pt2 - pts) + coefs[2 + up]*(ps2 + psm))/sqr(sHat()); double intr = 0.25*(coefs[4 + up]*pt2 + coefs[6 + up]*pts + coefs[8 + up]*ps2 + coefs[10 + up]*psm)* (sHat() - mZ2)/(sHat()*(sqr(sHat() - mZ2) + mZ2*GZ2)); lastBW = 0.25*(coefs[12 + up]*pt2 + coefs[14 + up]*pts + coefs[16 + up]*ps2 + coefs[18 + up]*psm)/ (sqr(sHat() - mZ2) + mZ2*GZ2); double alphaS = SM().alphaS(scale()); //int Nf = SM().Nf(scale()); int Nf = 3; DVector save; meInfo(save << lastCont << lastBW); return (lastCont + intr + lastBW)*sqr(SM().alphaEM(scale()))* (1.0 + alphaS/Constants::pi + (1.986-0.115*Nf)*sqr(alphaS/Constants::pi)); } Selector MEqq2gZ2ll::diagrams(const DiagramVector & diags) const { if ( lastXCombPtr() ) { lastCont = meInfo()[0]; lastBW = meInfo()[1]; } Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) { if ( diags[i]->id() == -1 ) sel.insert(lastCont, i); else if ( diags[i]->id() == -2 ) sel.insert(lastBW, i); } return sel; } Selector MEqq2gZ2ll::colourGeometries(tcDiagPtr diag) const { static ColourLines c("1 -2"); Selector sel; sel.insert(1.0, &c); return sel; } IBPtr MEqq2gZ2ll::clone() const { return new_ptr(*this); } IBPtr MEqq2gZ2ll::fullclone() const { return new_ptr(*this); } -void MEqq2gZ2ll::doinit() throw(InitException) { +void MEqq2gZ2ll::doinit() { double C = sqr(4.0*Constants::pi)/3.0; double SW2 = SM().sin2ThetaW(); double SW4 = sqr(SW2); double SW6 = SW2*SW4; double SW8 = SW2*SW6; double CW2 = 1.0 - SW2; coefs[0] = 16.0*C; coefs[1] = 64.0*C; coefs[2] = 8.0*C; coefs[3] = 32.0*C; C /= (CW2*SW2); coefs[4] = 4.0*(32.0*SW4 - 32.0*SW2 + 6.0)*C; coefs[5] = 8.0*(64.0*SW4 - 40.0*SW2 + 6.0)*C; coefs[6] = -4.0*(32.0*SW4 - 32.0*SW2 + 12.0)*C; coefs[7] = -8.0*(64.0*SW4 - 40.0*SW2 + 12.0)*C; coefs[8] = 4.0*(16.0*SW4 - 16.0*SW2 + 6.0)*C; coefs[9] = 8.0*(32.0*SW4 - 20.0*SW2 + 6.0)*C; coefs[10] = 4.0*(16.0*SW4 - 16.0*SW2 + 3.0)*C; coefs[11] = 8.0*(32.0*SW4 - 20.0*SW2 + 3.0)*C; C /= (CW2*SW2); coefs[12] = ( 64.0*SW8 - 128.0*SW6 + 128.0*SW4 - 48.0*SW2 + 9.0)*C; coefs[13] = (256.0*SW8 - 320.0*SW6 + 200.0*SW4 - 60.0*SW2 + 9.0)*C; coefs[14] = -( 64.0*SW8 - 128.0*SW6 + 176.0*SW4 - 96.0*SW2 + 18.0)*C; coefs[15] = -(256.0*SW8 - 320.0*SW6 + 296.0*SW4 - 120.0*SW2 + 18.0)*C; coefs[16] = ( 32.0*SW8 - 64.0*SW6 + 88.0*SW4 - 48.0*SW2 + 9.0)*C; coefs[17] = (128.0*SW8 - 160.0*SW6 + 148.0*SW4 - 60.0*SW2 + 9.0)*C; coefs[18] = ( 32.0*SW8 - 64.0*SW6 + 28.0*SW4 - 6.0*SW2)*C; coefs[19] = (128.0*SW8 - 160.0*SW6 + 64.0*SW4 - 12.0*SW2)*C; tcPDPtr Z0 = getParticleData(ParticleID::Z0); mZ2 = sqr(Z0->mass()); GZ2 = sqr(Z0->width()); ME2to2QCD::doinit(); } void MEqq2gZ2ll::persistentOutput(PersistentOStream & os) const { os << coefs << ounit(mZ2, GeV2) << ounit(GZ2, GeV2) << lastCont << lastBW; } void MEqq2gZ2ll::persistentInput(PersistentIStream & is, int) { is >> coefs >> iunit(mZ2, GeV2) >> iunit(GZ2, GeV2) >> lastCont >> lastBW; } ClassDescription MEqq2gZ2ll::initMEqq2gZ2ll; void MEqq2gZ2ll::Init() { static ClassDocumentation documentation ("The ThePEG::MEqq2gZ2ll class implements the full" "\\f$q\\bar{q}\\rightarrow\\gamma/Z^0\\rightarrow\\ell^+\\ell^-\\f$ " "matrix element including the interference terms."); } diff --git a/HerwigTest/MEqq2gZ2ll.h b/HerwigTest/MEqq2gZ2ll.h --- a/HerwigTest/MEqq2gZ2ll.h +++ b/HerwigTest/MEqq2gZ2ll.h @@ -1,128 +1,128 @@ // -*- C++ -*- #ifndef ThePEG_MEqq2gZ2ll_H #define ThePEG_MEqq2gZ2ll_H // // This is the declaration of the MEqq2gZ2ll class. // // CLASSDOC SUBSECTION Description: // // The MEqq2gZ2ll class implements the // e+e- -> gamma/Z0 -> // q+qbar matrix element. Both the continuum and // Z0 pole term as well as the interference term is // included. Although not a strict QCD matrix element the cass // inherits from ME2to2Base, mainly to inherit the // parameter for the number of active quark flavours. // // CLASSDOC SUBSECTION See also: // // ME2to2QCD.h. // #include "ThePEG/MatrixElement/ME2to2QCD.h" using namespace ThePEG; namespace Herwig { class MEqq2gZ2ll: public ME2to2QCD { public: MEqq2gZ2ll(); MEqq2gZ2ll(const MEqq2gZ2ll &); virtual ~MEqq2gZ2ll(); // Standard ctors and dtor public: virtual unsigned int orderInAlphaS() const; virtual unsigned int orderInAlphaEW() const; // Return the order in respective couplings in which this matrix // element is given. Returns 0 and 2 respectively. virtual double me2() const; // Return the matrix element for the kinematical configuation // previously provided by the last call to setKinematics(). virtual void getDiagrams() const; // Add all possible diagrams with the add() function. virtual Selector colourGeometries(tcDiagPtr diag) const; // Return a Selector with possible colour geometries for the selected // diagram weighted by their relative probabilities. virtual Selector diagrams(const DiagramVector &) const; // Weight the given diagrams with their relative probabilities. virtual Energy2 scale() const; // Return the scale associated with the last set phase space point. public: void persistentOutput(PersistentOStream &) const; void persistentInput(PersistentIStream &, int); // Standard functions for writing and reading from persistent streams. static void Init(); // Standard Init function used to initialize the interface. protected: virtual IBPtr clone() const; virtual IBPtr fullclone() const; // Standard clone methods - virtual void doinit() throw(InitException); + virtual void doinit(); // Standard Interfaced virtual functions. protected: vector coefs; // Constants for the different terms set from the StandardModel in // the init() function. Energy2 mZ2; Energy2 GZ2; // The mass squared and width squared of the Z0. mutable double lastCont; mutable double lastBW; // The last continuum and Breit-Wigner terms to be used to select // primary diagram. private: static ClassDescription initMEqq2gZ2ll; MEqq2gZ2ll & operator=(const MEqq2gZ2ll &); // Private and non-existent assignment operator. }; } namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ template <> struct BaseClassTrait: public ClassTraitsType { typedef ME2to2QCD NthBase; }; template <> struct ClassTraits : public ClassTraitsBase { static string className() { return "/Herwig++/MEqq2gZ2ll"; } static string library() { return "MEqq2gZ2ll.so"; } }; /** @endcond */ } #include "MEqq2gZ2ll.icc" #endif /* ThePEG_MEqq2gZ2ll_H */