diff --git a/src/dyturbo.C b/src/dyturbo.C index 913d2bf..807c066 100644 --- a/src/dyturbo.C +++ b/src/dyturbo.C @@ -1,385 +1,434 @@ #ifndef dyturbo_C #define dyturbo_C /** * @file dyturbo.C * Steering class for Drell-Yan calculation. * * @brief A brief description * * @author Jakub Cuth * @date 2016-09-17 */ #include "dyturbo.h" #include "cubacall.h" #include "clock_real.h" #include "settings.h" #include "interface.h" #include "ctintegr.h" #include "phasespace.h" #include "rapint.h" #include "HistoHandler.h" using DYTurbo::PrintTable::Col3; using DYTurbo::PrintTable::Col4; bool DYTurbo::HasOnlyVegas = false; bool DYTurbo::isDryRun = false; -namespace DYTurbo { - - // definition of data member - Term subtotal; - TermList ActiveTerms; - BoundariesList ActiveBoundaries; - BoundIterator last_bounds; - - // Term - - void Term::last_reset() { - last_int.assign(opts.totpdf,0.); - last_err2=0; - last_time = clock_real(); - } - - void Term::RunIntegration(){ - double err; - last_reset(); - /// @todo specialized (need to reincorporate) - if ( integrate == resintegr2d ){ - if (opts.resumcpp) rapint::cache(phasespace::ymin, phasespace::ymax); - else cacheyrapint_(phasespace::ymin, phasespace::ymax); - } - // run - if (isDryRun){ - // this is for testing interface - last_int[0] = 1; err = 1; - } else { - integrate(last_int,err); - } - // - last_time = clock_real()-last_time; - last_err2 += err*err; - // save results to histograms - if (!isVegas){ - for (size_t ivar = 0; ivar < last_int.size(); ++ivar) { - HistoHandler::SetVariation(ivar); - HistoHandler::FillResult(last_int[ivar],sqrt(last_err2)); - } - } - // cumulate - total_time+=last_time; - total_int+=last_int[0]; - total_err2+=last_err2; - // cumulate integral to subtotal - subtotal.last_int[0] += last_int[0]; - subtotal.total_int += last_int[0]; - // cumulate error to subtotal - subtotal.last_err2 += last_err2; - subtotal.total_err2 += last_err2; - } - - void Term::Print(){ - printf("%24s:%s",name.c_str(), description.c_str()); - } - - // Term iterator - - TermIterator::TermIterator() : icurrent(0) { } - - bool TermIterator::IsEnd(){ - return icurrent==ActiveTerms.size(); - } - - TermIterator & TermIterator::operator++(){ - icurrent++; - return (*this); - } - - Term & TermIterator::operator*(){ - return ActiveTerms[icurrent]; - } - - // Boundary iterator - - BoundIterator::BoundIterator(){ - // set first boundary - isFirst=true; - current.clear(); - if (!ActiveBoundaries.empty()) - for (size_t i = 0; i < N_boundaries; ++i) { - current.push_back(ActiveBoundaries[i].begin()); - } +namespace DYTurbo +{ + // definition of data member + Term subtotal; + TermList ActiveTerms; + BoundariesList ActiveBoundaries; + BoundIterator last_bounds; + + //Term + + void Term::last_reset() + { + last_int.assign(opts.totpdf,0.); + last_err2=0; + last_time = clock_real(); + } + + void Term::RunIntegration() + { + double err; + last_reset(); + /// @todo specialized (need to reincorporate) + if ( integrate == resintegr2d ){ + if (opts.resumcpp) rapint::cache(phasespace::ymin, phasespace::ymax); + else cacheyrapint_(phasespace::ymin, phasespace::ymax); } - - bool BoundIterator::IsEnd(){ - if (isFirst) return false; - for (size_t i = 0; i < N_boundaries; ++i) - if (current[i] != ActiveBoundaries[i].begin() ) - return false; - return true; + // run + if (isDryRun) // this is for testing interface + { + last_int[0] = 1; err = 1; + } + else + integrate(last_int,err); + // + if (opts.ptbinwidth) {last_int[0] /= (phasespace::qtmax-phasespace::qtmin); err /= (phasespace::qtmax-phasespace::qtmin);} + if (opts.ybinwidth) {last_int[0] /= (phasespace::ymax-phasespace::ymin); err /= (phasespace::ymax-phasespace::ymin);} + if (opts.mbinwidth) {last_int[0] /= (phasespace::mmax-phasespace::mmin); err /= (phasespace::mmax-phasespace::mmin);} + last_time = clock_real()-last_time; + last_err2 += err*err; + // save results to histograms + if (!isVegas) + { + if (opts.makehistos) + for (size_t ivar = 0; ivar < last_int.size(); ++ivar) + { + HistoHandler::SetVariation(ivar); + HistoHandler::FillResult(last_int[ivar],sqrt(last_err2)); + } + } + // cumulate + total_time+=last_time; + total_int+=last_int[0]; + total_err2+=last_err2; + // cumulate integral to subtotal + subtotal.last_int[0] += last_int[0]; + subtotal.total_int += last_int[0]; + // cumulate error to subtotal + subtotal.last_err2 += last_err2; + subtotal.total_err2 += last_err2; + + //// cumulate integral to bintotal + //bintotal.last_int[0] += bin_int[0]; + //bintotal.total_int += last_int[0]; + // cumulate error to subtotal + //bintotal.last_err2 += last_err2; + //bintotal.total_err2 += last_err2; + } + + void Term::Print() + { + printf("%24s:%s",name.c_str(), description.c_str()); + } + + // Term iterator + TermIterator::TermIterator() : icurrent(0) { } + + bool TermIterator::IsEnd() + { + return icurrent==ActiveTerms.size(); + } + + TermIterator & TermIterator::operator++() + { + icurrent++; + return (*this); + } + + Term & TermIterator::operator*() + { + return ActiveTerms[icurrent]; + } + + // Boundary iterator + BoundIterator::BoundIterator() + { + // set first boundary + isFirst=true; + current.clear(); + if (!ActiveBoundaries.empty()) + for (size_t i = 0; i < N_boundaries; ++i) + current.push_back(ActiveBoundaries[i].begin()); + + } + + bool BoundIterator::IsEnd() + { + if (isFirst) return false; + for (size_t i = 0; i < N_boundaries; ++i) + if (current[i] != ActiveBoundaries[i].begin() ) + return false; + return true; + } + + BoundIterator & BoundIterator::operator++() + { + for (int i = N_boundaries-1; i >= 0; --i) { // start from last + current[i]++; // increase + // check if it's equal to previous to last (need two numbers as boundaries) + if(current[i]!=ActiveBoundaries[i].end()-1) break; // go to return + else current[i]=ActiveBoundaries[i].begin(); // is previous to last } - - BoundIterator & BoundIterator::operator++(){ - for (int i = N_boundaries-1; i >= 0; --i) { // start from last - current[i]++; // increase - // check if it's equal to previous to last (need two numbers as boundaries) - if(current[i]!=ActiveBoundaries[i].end()-1) break; // go to return - else current[i]=ActiveBoundaries[i].begin(); // is previous to last - } - /** - * @attention After looping through all of boundary bins it points back - * to beginning. Therefore it is necessary to set `isFirst=false` - */ - if (isFirst) isFirst=false; - return (*this); - } - - void BoundIterator::Print(){ - printf ("Boundaries "); - printf ( "| M : %f -- %f " , *current [ b_M ] , * ( current [ b_M ] +1 ) ) ; - printf ( "| Qt : %f -- %f " , *current [ b_QT ] , * ( current [ b_QT ] +1 ) ) ; - printf ( "| Y : %f -- %f " , *current [ b_Y ] , * ( current [ b_Y ] +1 ) ) ; - printf ( "| CsTh : %f -- %f " , *current [ b_CsTh ] , * ( current [ b_CsTh ] +1 ) ) ; - printf ( "\n"); - } - - - - //! Internal flag only for test purposes - bool TestAllTerms=false; - /** - * @brief Add term to active terms if is requested. - * - * Also it is checked if we are only using Vegas for integration. TYhi - * - * @param isActive If this false then skip all stuff. - * @param fun The pointer to integration function, which should be called for calculation. - * @param name Set the name of term - * @param is_vegas For checking if we need Integration mode. - * - * @return It returns newly added term so we can define description by using stream operator. + * @attention After looping through all of boundary bins it points back + * to beginning. Therefore it is necessary to set `isFirst=false` */ - Term & AddTermIfActive(const bool &isActive, void (* fun)(VecDbl &val,double &err), const String &name, const bool &is_vegas ){ - subtotal.description = ""; - if (!isActive && !TestAllTerms) return subtotal; - ActiveTerms.push_back(Term()); - ActiveTerms.back().name = name; - ActiveTerms.back().description = ""; - ActiveTerms.back().integrate = fun; - ActiveTerms.back().isVegas = is_vegas; - HasOnlyVegas&=is_vegas; - return ActiveTerms.back(); - } + if (isFirst) isFirst=false; + return (*this); + } + + void BoundIterator::Print(){ + printf ("Boundaries "); + printf ( "| M : %f -- %f " , *current [ b_M ] , * ( current [ b_M ] +1 ) ) ; + printf ( "| Qt : %f -- %f " , *current [ b_QT ] , * ( current [ b_QT ] +1 ) ) ; + printf ( "| Y : %f -- %f " , *current [ b_Y ] , * ( current [ b_Y ] +1 ) ) ; + printf ( "| CsTh : %f -- %f " , *current [ b_CsTh ] , * ( current [ b_CsTh ] +1 ) ) ; + printf ( "\n"); + } + + //! Internal flag only for test purposes + bool TestAllTerms=false; + + /** + * @brief Add term to active terms if is requested. + * + * Also it is checked if we are only using Vegas for integration. TYhi + * + * @param isActive If this false then skip all stuff. + * @param fun The pointer to integration function, which should be called for calculation. + * @param name Set the name of term + * @param is_vegas For checking if we need Integration mode. + * + * @return It returns newly added term so we can define description by using stream operator. + */ + Term & AddTermIfActive(const bool &isActive, void (* fun)(VecDbl &val,double &err), const String &name, const bool &is_vegas ){ + subtotal.description = ""; + if (!isActive && !TestAllTerms) return subtotal; + ActiveTerms.push_back(Term()); + ActiveTerms.back().name = name; + ActiveTerms.back().description = ""; + ActiveTerms.back().integrate = fun; + ActiveTerms.back().isVegas = is_vegas; + HasOnlyVegas&=is_vegas; + return ActiveTerms.back(); + } template Term & Term::operator<<(const Streamable &data){ SStream strm; strm << data; description += strm.str(); return (*this); } void AddTerms(){ ActiveTerms.clear(); HasOnlyVegas=true; const bool isVegas=true; const bool isNotVegas=false; int w0=25; int w1=30; int w2=12; int w3=12; string name; // finite born bool fixed_born = opts.doBORN && opts.fixedorder; if (fixed_born || TestAllTerms) { name="Fixed born"; AddTermIfActive ( opts.bornint2d , bornintegr2d , name, isNotVegas ) << Col3( "cuhre (dm, dpt)" , "iter =" , opts.niterBORN ); AddTermIfActive ( opts.bornintvegas4d , bornintegrMC4d , name, isVegas ) << Col3( "vegas 4D" , "ncalls =" , opts.vegasncallsBORN ); AddTermIfActive ( opts.bornintvegas6d , bornintegrMC6d , name, isVegas ) << Col3( "vegas 6D" , "ncalls =" , opts.vegasncallsBORN ); } // resummation bool resum_born = opts.doBORN && !opts.fixedorder; if (resum_born || TestAllTerms) { name="Resummation"; AddTermIfActive ( opts.resint1d , resintegr1d , name , isNotVegas ) << Col3 ( "cuhre (dpt)" , "iter =" , opts.niterBORN ) << Col4 ( "","gauss (dy)" , "nodes =" , opts.yrule ) << Col4 ( "","" , "intervals =" , opts.yintervals ) << Col3 ( "","analytical (dpt)","" ); AddTermIfActive ( opts.resint2d , resintegr2d , name , isNotVegas ) << Col3 ( "cuhre (dm, dpt)" , "iter =" , opts.niterBORN ) << Col4 ( "","gauss (dy)" , "nodes =" , opts.yrule ) << Col4 ( "","" , "intervals =" , opts.yintervals ); AddTermIfActive ( opts.resint3d , resintegr3d , name , isNotVegas ) << Col3 ( "cuhre (dm, dpt, dy)" , "iter =" , opts.niterBORN ); AddTermIfActive ( opts.resintvegas , resintegrMC , name , isVegas ) << Col3 ( "vegas" , "ncalls =" , opts.vegasncallsBORN ); } // CT if (opts.doCT || TestAllTerms) { name="Counter term"; AddTermIfActive ( opts.ctint1d , ctintegr1d , name , isNotVegas ) << Col3 ( "cuhre (dm)" , "iter =" , opts.niterCT ) << Col4 ( "","gauss (dpt)" , "nodes =" , opts.qtrule ) << Col4 ( "","" , "intervals =" , opts.qtintervals ); AddTermIfActive ( opts.ctint2d , ctintegr2d , name , isNotVegas ) << Col3 ( "cuhre (dm, dy)" , "iter =" , opts.niterCT ) << Col4 ( "","gauss (dpt)" , "nodes =" , 20 ) << Col4 ( "","" , "intervals =" , 5 ); AddTermIfActive ( opts.ctint3d , ctintegr3d , name , isNotVegas ) << Col3 ( "cuhre (dm, dpt, dy)" , "iter =" , opts.niterCT ); AddTermIfActive ( opts.ctintvegas6d , ctintegrMC , name , isVegas ) << Col3 ( "vegas 6D" , "ncalls =" , opts.vegasncallsCT ); AddTermIfActive ( opts.ctintvegas8d , ctintegr , name , isVegas ) << Col3 ( "vegas 8D" , "ncalls =" , opts.vegasncallsCT ); } // VJ finite bool vj_finite = opts.doVJ; name="V+J"; if (vj_finite || TestAllTerms){ AddTermIfActive ( opts.vjint3d , vjintegr3d , name , isNotVegas ) << Col3 ( "cuhre (dm, dpt, dy)" , "iter =" , opts.niterVJ ); // VJ LO if (opts.order == 1) { AddTermIfActive ( opts.vjint5d && opts.order == 1 , vjlointegr5d , name , isNotVegas ) << Col3 ( "cuhre 5D???" , "iter =" , opts.niterVJ ); AddTermIfActive ( opts.vjintvegas7d && opts.order == 1 , vjlointegr7d , name , isVegas ) << Col3 ( "vegas 7D" , "ncalls =" , opts.vegasncallsVJLO ); //phase space improved MCFM integration //AddTermIfActive ( opts.vjintvegas7d && opts.order == 1 , vjlointegr , name , isVegas ) << Col3 ( "vegas 7D" , "ncalls =" , opts.vegasncallsVJLO ); //original MCFM integration } // VJ NLO if (!opts.vjint3d && opts.order == 2) { AddTermIfActive ( opts.doVJREAL , vjrealintegr , "V+J Real" , isVegas) << Col3 ( "vegas" , "ncalls =" , opts.vegasncallsVJREAL ); AddTermIfActive ( opts.doVJVIRT , vjvirtintegr , "V+J Virtual" , isVegas) << Col3 ( "vegas" , "ncalls =" , opts.vegasncallsVJVIRT ); } } } /** * @brief Check for Integration mode and add boundary. * * If we have only Vegas calculations we can run with simple boundaries. * This means we will only pick up one bin with lowest and hightest value * as boundaries. * * @param ib Index of variable (\ref BoundaryIndex can be used as input) * @param name Name of boundary. Will be used for printing. * @param bins All boundaries requested from user. * @return Description of return. */ void AddBoundary(size_t ib, string name, VecDbl &bins ){ ActiveBoundaries[ib].name = name; ActiveBoundaries[ib].data.clear(); if (HasOnlyVegas && !opts.force_binsampling){ // simple boundary mode ActiveBoundaries[ib].data.push_back(bins.front()); ActiveBoundaries[ib].data.push_back(bins.back()); } else { // integration mode ActiveBoundaries[ib].data = bins; } } /** * @brief Adding all boundaries * * This have to be done after adding all terms (becaus of check for Vegas * only). * * If adding new boundary definition please also add new enum item into \ref BoundaryIndex. * * @param _inArg Description of param * @return Description of return. */ - void AddBoundaries(){ - ActiveBoundaries.resize(N_boundaries); - VecDbl costh{opts.costhmin , opts.costhmax}; - AddBoundary ( b_M , "q2" , bins.mbins ) ; - AddBoundary ( b_Y , "y" , bins.ybins ) ; - AddBoundary ( b_QT , "qT" , bins.qtbins ) ; - AddBoundary ( b_CsTh , "costh" , costh ); + void AddBoundaries() + { + ActiveBoundaries.resize(N_boundaries); + VecDbl costh{opts.costhmin , opts.costhmax}; + AddBoundary ( b_M , "q2" , bins.mbins ) ; + AddBoundary ( b_Y , "y" , bins.ybins ) ; + AddBoundary ( b_QT , "qT" , bins.qtbins ) ; + AddBoundary ( b_CsTh , "costh" , costh ); } /** * @brief First run of resummation and counter term. * * If using the DYRES approximation for PDFs, make sure that the PDF fit * is initialised in the same way Need to throw a random point according * to a breit wigner, which is used to determine xtauf in the PDF fit * */ void WarmUpResummation(){ double costh, m, qt, y; int mode = 0; if (opts.doBORN || opts.doCT){ if (opts_.approxpdf_ == 1) { srand(opts.rseed); double wsqmin = pow(phasespace::mmin ,2); double wsqmax = pow(phasespace::mmax ,2); double x1=((double)rand()/(double)RAND_MAX); double m2,wt; breitw_(x1,wsqmin,wsqmax,opts.rmass,opts.rwidth,m2,wt); printf( "Initialise PDF fit with mass = %f xtauf = %f \n", sqrt(m2), m2/opts.sroot); costh = 0.; m = sqrt(m2); qt = 1; y = 0; for (int ipdf=0; ipdf &val, vector &err) + { + //Set up integration terms and bins boundaries + DYTurbo::WarmUp(); + + //Loop on phase space bins + for ( DYTurbo::BoundIterator bounds; !bounds.IsEnd(); ++bounds) + { + DYTurbo::SetBounds(bounds); + + //Loop on active terms + for (DYTurbo::TermIterator term; !term.IsEnd(); ++term) + { + (*term).RunIntegration(); + } + + val.push_back(subtotal.last_int[0]); + err.push_back(sqrt(subtotal.last_err2)); + //cout << "dyturbo " << val.size() << " " << subtotal.last_int[0] << " " << sqrt(subtotal.last_err2) << endl; + subtotal.last_reset(); + } + DYTurbo::Terminate(); + } }; #endif /* ifndef dyturbo_C */ diff --git a/src/dyturbo.h b/src/dyturbo.h index f98dc1e..2620c8c 100644 --- a/src/dyturbo.h +++ b/src/dyturbo.h @@ -1,319 +1,328 @@ #ifndef dyturbo_H #define dyturbo_H /** * @file dyturbo.h * Steering class for Drell-Yan calculation. * * @brief A brief description * * @author Jakub Cuth * @date 2016-09-17 */ #include using std::setw; using std::setprecision; #include "handy_typdefs.h" /** * @brief Interface to properly control calculation. * */ -namespace DYTurbo { +namespace DYTurbo +{ /** * @brief Flag for if we need Integration mode. * * If we have only Vegas integration we don't need to run per each bin and * boundaries can be simplified. Only one bin from lowest to highest * value. */ - extern bool HasOnlyVegas; - - //! Debug flag for testing code. - extern bool isDryRun; - - - //forward - struct Term; - struct Boundaries; - struct BoundIterator; - - /** @defgroup MainFunctions Main DYTurbo functions. - * @{ - */ + extern bool HasOnlyVegas; + + //! Debug flag for testing code. + extern bool isDryRun; + + + //forward + struct Term; + struct Boundaries; + struct BoundIterator; + + /** @defgroup MainFunctions Main DYTurbo functions. + * @{ + */ + + /** + * @brief Initialization of program. + * + * It parse command line and input file and set interfaced values. It's + * also calling initialization of submodules. + */ + void Init(int argc, char * argv[]); + + //Initialise constants + void init_const(); + /** + * @brief Initialization of submodules. + */ + void init_params(); + + //release memory allocated by init_params() + void release(); + + /** + * @brief Warming up integration. + * + * Some of terms needed to be pre-run before actual integration. This done + * inside this function. It also prints out current integration settings. + */ + void WarmUp(); + + /** + * @brief Set integration bounds. + * + * The values from Boundaries are read and provided to \ref phasespace. + * + * @attention This needs to be called before calling the \ref RunIntegration. + */ + void SetBounds(BoundIterator); + + /** + * @brief Clean up and save results. + */ + void Terminate(); + + void compute(std::vector &val, std::vector &err); + + namespace PrintTable { + void Init(); + void Settings(); + void Header() ; + void Bounds(bool use_full_bound=false) ; + void Result(const Term &term,bool printGrandTotal=false) ; + void ResultSubTotal(bool is_grandtotal=false) ; + void ResultGrandTotal() ; + void Footer() ; + } + /** @} */ + + /** + * @brief Class for keeping information about active term. + * + * This class contains name, description and function to be called for + * selected term. Terms are then stored in TermList and looped in main + * program. + */ + struct Term { + String name = "dummy"; //!< Name of term will be used in header of table. + String description = ""; //!< Description is used while calling \ref Print function. + bool isVegas = true; //!< If term is not Vegas fill result. /** - * @brief Initialization of program. + * @brief Pointer to integration function. This par * - * It parse command line and input file and set interfaced values. It's - * also calling initialization of submodules. + * This will be called within \ref RunIntegration function. */ - void Init(int argc, char * argv[]); - /** - * @brief Initialization of submodules. - */ - void init_params(); + void (*integrate)(VecDbl &val,double &err) = NULL; - /** - * @brief Warming up integration. - * - * Some of terms needed to be pre-run before actual integration. This done - * inside this function. It also prints out current integration settings. - */ - void WarmUp(); + double total_int = 0; //!< Total cross-section summed from all bins. + double total_err2 = 0; //!< Squared uncertainties summed from all bins. + double total_time = 0; //!< Total time spend for all bins. + VecDbl last_int = {}; //!< Here is last cross-section stored for all variations (PDF). + double last_err2 =0; //!< Uncertainty of last integration. + double last_time =0; //!< Time spend in last iteration. /** - * @brief Set integration bounds. - * - * The values from Boundaries are read and provided to \ref phasespace. + * @brief Main functionality of this class. * - * @attention This needs to be called before calling the \ref RunIntegration. + * This will call integration, measure time for processing, store and + * cumulate the result. */ - void SetBounds(BoundIterator); + void RunIntegration(); /** - * @brief Clean up and save results. - */ - void Terminate(); - - namespace PrintTable { - void Init(); - void Settings(); - void Header() ; - void Bounds(bool use_full_bound=false) ; - void Result(const Term &term,bool printGrandTotal=false) ; - void ResultSubTotal(bool is_grandtotal=false) ; - void ResultGrandTotal() ; - void Footer() ; - } - /** @} */ - - /** - * @brief Class for keeping information about active term. + * @brief Print settings of this term to screen. * - * This class contains name, description and function to be called for - * selected term. Terms are then stored in TermList and looped in main - * program. + * This is used while printing \ref IntegrationSettings and format + * follow 4 column mode. */ - struct Term { - String name = "dummy"; //!< Name of term will be used in header of table. - String description = ""; //!< Description is used while calling \ref Print function. - bool isVegas = true; //!< If term is not Vegas fill result. - - /** - * @brief Pointer to integration function. This par - * - * This will be called within \ref RunIntegration function. - */ - void (*integrate)(VecDbl &val,double &err) = NULL; - - double total_int = 0; //!< Total cross-section summed from all bins. - double total_err2 = 0; //!< Squared uncertainties summed from all bins. - double total_time = 0; //!< Total time spend for all bins. - VecDbl last_int = {}; //!< Here is last cross-section stored for all variations (PDF). - double last_err2 =0; //!< Uncertainty of last integration. - double last_time =0; //!< Time spend in last iteration. - - /** - * @brief Main functionality of this class. - * - * This will call integration, measure time for processing, store and - * cumulate the result. - */ - void RunIntegration(); - - /** - * @brief Print settings of this term to screen. - * - * This is used while printing \ref IntegrationSettings and format - * follow 4 column mode. - */ - void Print(); - - //! Overloading stream operator to set description. - template Term & operator<<(const Streamable &data); - - //! Reset to zero values of last_int, last_err2 and last_time - void last_reset(); + void Print(); + + //! Overloading stream operator to set description. + template Term & operator<<(const Streamable &data); + + //! Reset to zero values of last_int, last_err2 and last_time + void last_reset(); + }; + + //! Helper instance for summing results from all terms. + extern Term subtotal; + + //! Container class for Terms + typedef std::vector TermList; + + /** + * @brief Helper instance storing currenlty active terms. + * + * This is where DYTurbo stores active terms. They are looped by \ref + * TermIterator. + */ + extern TermList ActiveTerms; + + /** + * @brief Helper struct for looping active terms. + * + * Term Iterator always loops instance \ref ActiveTerms from beginning. + */ + struct TermIterator { + TermIterator(); //!< By construction always starts at begining of ActiveTerms. + size_t icurrent; //!< Storing position in vector. + bool IsEnd(); //!< Test iterator is in the end. + TermIterator& operator++(); //! Increase operator simple increases \ref icurrent variable. + Term & operator*(); //!< Dereferencing operator uses `vector::operator[]` + }; + + + + //! Container of boundary values for one variable. + typedef VecDbl BoundaryValues; + //! Iterator of boundary values for one variable. + typedef BoundaryValues::iterator BoundValueItr; + + /** + * @brief Holder of boundaries for one variable + */ + struct Boundaries{ + std::string name =""; //!< Name of variable for printing purpose. + BoundaryValues data = {}; //!< Actual values of boundaries for this variable. + + inline BoundValueItr begin() { return data.begin();}; //!< Iterator interface to data. + inline BoundValueItr end() { return data.end(); }; //!< Iterator interface to data. + inline double front() { return data.front();}; //!< Iterator interface to data. + inline double back() { return data.back(); }; //!< Iterator interface to data. + inline size_t size() { return data.size(); }; //!< Iterator interface to data. + }; + + /** + * @brief Define position of variable inside \ref ActiveBoundaries + * + * When adding new variable please also add filling in \ref AddBoundaries + * and interface to phasespace boundaries in \ref SetBounds. + * + * @note Trick is to define one extra enum item in the end for counting number of + * defined items. + * + * @attention The `N_boundaries` must be always the last item. + * + * @todo b_Phi, b_LepPt + */ + enum BoundaryIndex + { + b_M=0, //!< Vector boson invariant mass + b_Y, //!< Vector boson rapidity + b_QT, //!< Vector boson transverse momentum + b_CsTh, //!< Cosine of longitudinal angle in Collins-Soper frame + N_boundaries //! This here to count number of boundaries. }; - //! Helper instance for summing results from all terms. - extern Term subtotal; - - //! Container class for Terms - typedef std::vector TermList; - - /** - * @brief Helper instance storing currenlty active terms. - * - * This is where DYTurbo stores active terms. They are looped by \ref - * TermIterator. - */ - extern TermList ActiveTerms; + /** + * @brief Container of all variable boundaries. + * + * Position of variable is set by BoundaryIndex. This is filled inside function \ref AddBoundaries. + */ + typedef std::vector BoundariesList; + + /** + * @brief Helper instance containing active boundaries. + * + * This is looped by BoundIterator. + */ + extern BoundariesList ActiveBoundaries; + + /** + * @brief Container of boundary value iterators + */ + typedef std::vector BoundValIterList; + + + /** + * @brief Iterating over all variable boundaries (bins) for integration. + */ + struct BoundIterator { + BoundIterator(); //!< After construction always point to first boundary. + bool IsEnd(); //!< Check we already looped through all boundaries. /** - * @brief Helper struct for looping active terms. + * @brief Go to next kinematic boundary. * - * Term Iterator always loops instance \ref ActiveTerms from beginning. - */ - struct TermIterator { - TermIterator(); //!< By construction always starts at begining of ActiveTerms. - size_t icurrent; //!< Storing position in vector. - bool IsEnd(); //!< Test iterator is in the end. - TermIterator& operator++(); //! Increase operator simple increases \ref icurrent variable. - Term & operator*(); //!< Dereferencing operator uses `vector::operator[]` - }; - - - - //! Container of boundary values for one variable. - typedef VecDbl BoundaryValues; - //! Iterator of boundary values for one variable. - typedef BoundaryValues::iterator BoundValueItr; - - /** - * @brief Holder of boundaries for one variable + * This function assures that all combinations of boundaries are taken. + * For more detail see implementation. */ - struct Boundaries{ - std::string name =""; //!< Name of variable for printing purpose. - BoundaryValues data = {}; //!< Actual values of boundaries for this variable. - - inline BoundValueItr begin() { return data.begin();}; //!< Iterator interface to data. - inline BoundValueItr end() { return data.end(); }; //!< Iterator interface to data. - inline double front() { return data.front();}; //!< Iterator interface to data. - inline double back() { return data.back(); }; //!< Iterator interface to data. - inline size_t size() { return data.size(); }; //!< Iterator interface to data. - }; + BoundIterator& operator++(); - /** - * @brief Define position of variable inside \ref ActiveBoundaries - * - * When adding new variable please also add filling in \ref AddBoundaries - * and interface to phasespace boundaries in \ref SetBounds. - * - * @note Trick is to define one extra enum item in the end for counting number of - * defined items. - * - * @attention The `N_boundaries` must be always the last item. - * - * @todo b_Phi, b_LepPt - */ - enum BoundaryIndex { - b_M=0, //!< Vector boson invariant mass - b_Y, //!< Vector boson rapidity - b_QT, //!< Vector boson transverse momentum - b_CsTh, //!< Cosine of longitudinal angle in Collins-Soper frame - N_boundaries //! This here to count number of boundaries. - }; + //! Print to screen. Used by `PrintTable::Boundary`. + void Print(); /** - * @brief Container of all variable boundaries. - * - * Position of variable is set by BoundaryIndex. This is filled inside function \ref AddBoundaries. + * @brief Current lower boundary. + * @param ib Index of variable (\ref BoundaryIndex can be used as input) + * @return Lower boundary of selected variable. */ - typedef std::vector BoundariesList; + inline double loBound(size_t ib){return *current[ib] ;} /** - * @brief Helper instance containing active boundaries. - * - * This is looped by BoundIterator. + * @brief Current upper boundary. + * @param ib Index of variable (\ref BoundaryIndex can be used as input) + * @return Upper boundary of selected variable. */ - extern BoundariesList ActiveBoundaries; + inline double hiBound(size_t ib){return *(current[ib]+1) ;} /** - * @brief Container of boundary value iterators + * @brief List of pointers to all variable boundary values. */ - typedef std::vector BoundValIterList; - - + BoundValIterList current; + /** - * @brief Iterating over all variable boundaries (bins) for integration. + * @brief Control flag. + * + * Set to true in beginning and after first increment is changed to + * false. */ - struct BoundIterator { - BoundIterator(); //!< After construction always point to first boundary. - bool IsEnd(); //!< Check we already looped through all boundaries. - - /** - * @brief Go to next kinematic boundary. - * - * This function assures that all combinations of boundaries are taken. - * For more detail see implementation. - */ - BoundIterator& operator++(); - - //! Print to screen. Used by `PrintTable::Boundary`. - void Print(); - - /** - * @brief Current lower boundary. - * @param ib Index of variable (\ref BoundaryIndex can be used as input) - * @return Lower boundary of selected variable. - */ - inline double loBound(size_t ib){return *current[ib] ;} - - /** - * @brief Current upper boundary. - * @param ib Index of variable (\ref BoundaryIndex can be used as input) - * @return Upper boundary of selected variable. - */ - inline double hiBound(size_t ib){return *(current[ib]+1) ;} - - /** - * @brief List of pointers to all variable boundary values. - */ - BoundValIterList current; - - /** - * @brief Control flag. - * - * Set to true in beginning and after first increment is changed to - * false. - */ - bool isFirst; - }; + bool isFirst; + }; - //! Helper bounds: remember last bounds for printing. - extern BoundIterator last_bounds; - - - namespace PrintTable{ - //! Helper structs for printing: Four column - const int colw1 = 25; - const int colw2 = 20; - const int colw3 = 15; - const int colw4 = 12; - - struct Col4 { - String data; - template Col4( S1 col1, S2 col2, S3 col3, S4 col4, int prec=6){ - SStream tmp; - tmp << setprecision(prec) << setw(colw1) << col1; - tmp << setprecision(prec) << setw(colw2) << col2; - tmp << setprecision(prec) << setw(colw3) << col3; - tmp << setprecision(prec) << setw(colw4) << col4; - tmp << '\n'; - data = tmp.str(); - } - }; - inline std::ostream & operator<< (std::ostream & strm, const Col4 &col){ strm << col.data; return strm; } - - //! Helper structs for printing: Four column, skip first - struct Col3{ - String data; - template Col3(S2 col2, S3 col3, S4 col4){ - SStream tmp; - tmp << setw(colw2) << col2; - tmp << setw(colw3) << col3; - tmp << setw(colw4) << col4; - tmp << '\n'; - data = tmp.str(); - } - }; - inline std::ostream & operator<< (std::ostream & strm, const Col3 &col){ strm << col.data; return strm; } - } - + //! Helper bounds: remember last bounds for printing. + extern BoundIterator last_bounds; + + + namespace PrintTable{ + //! Helper structs for printing: Four column + const int colw1 = 25; + const int colw2 = 20; + const int colw3 = 15; + const int colw4 = 12; + + struct Col4 { + String data; + template Col4( S1 col1, S2 col2, S3 col3, S4 col4, int prec=6){ + SStream tmp; + tmp << setprecision(prec) << setw(colw1) << col1; + tmp << setprecision(prec) << setw(colw2) << col2; + tmp << setprecision(prec) << setw(colw3) << col3; + tmp << setprecision(prec) << setw(colw4) << col4; + tmp << '\n'; + data = tmp.str(); + } + }; + inline std::ostream & operator<< (std::ostream & strm, const Col4 &col){ strm << col.data; return strm; } + + //! Helper structs for printing: Four column, skip first + struct Col3{ + String data; + template Col3(S2 col2, S3 col3, S4 col4){ + SStream tmp; + tmp << setw(colw2) << col2; + tmp << setw(colw3) << col3; + tmp << setw(colw4) << col4; + tmp << '\n'; + data = tmp.str(); + } + }; + inline std::ostream & operator<< (std::ostream & strm, const Col3 &col){ strm << col.data; return strm; } + } } #endif /* ifndef dyturbo_H */