diff --git a/MatrixElement/FxFx/FxFxAnalysis.cc b/MatrixElement/FxFx/FxFxAnalysis.cc deleted file mode 100644 --- a/MatrixElement/FxFx/FxFxAnalysis.cc +++ /dev/null @@ -1,389 +0,0 @@ -// -*- C++ -*- -// -// This is the implementation of the non-inlined, non-templated member -// functions of the FxFxAnalysis class. -// - -#include "FxFxAnalysis.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/ParVector.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Interface/Switch.h" -#include "ThePEG/Persistency/PersistentOStream.h" -#include "ThePEG/Persistency/PersistentIStream.h" -#include "ThePEG/Vectors/HepMCConverter.h" -#include "ThePEG/Config/HepMCHelper.h" -#include "ThePEG/EventRecord/Event.h" -#include "ThePEG/EventRecord/SubProcess.h" -#include "ThePEG/EventRecord/SubProcessGroup.h" -#include "ThePEG/Repository/EventGenerator.h" -#include "ThePEG/Repository/CurrentGenerator.h" -#include "HepMC/GenEvent.h" -#include "Rivet/AnalysisHandler.hh" -#include "Rivet/Tools/Logging.hh" -#include "HepMC/IO_AsciiParticles.h" -#include "HepMC/IO_GenEvent.h" -#include "ThePEG/Utilities/XSecStat.h" - -using namespace ThePEG; - -FxFxAnalysis::FxFxAnalysis() - : _remnantId(82), _format(1),_unitchoice(), - _geneventPrecision(16), debug(false),useoptweights(false), normoptweights(false), _rivet(), _nevent(0) {} - -HepMC::GenEvent * FxFxAnalysis::makeEvent(tEventPtr event, tSubProPtr sub, long no, - Energy eUnit, Length lUnit, - CrossSection xsec, CrossSection xsecErr) const { - //convert the event from the Herwig format to the HepMC format and write it to the common block - HepMC::GenEvent * ev = HepMCConverter::convert(*event, false,eUnit, lUnit); - - //reset the event - HepMCTraits::resetEvent(ev, no, event->weight()*sub->groupWeight(), event->optionalWeights()); - - //set the cross section - HepMCTraits::setCrossSection(*ev,xsec/picobarn, - xsecErr/picobarn); - return ev; -} - -HepMC::GenEvent * FxFxAnalysis::makeEventW(tEventPtr event, long no, - Energy eUnit, Length lUnit, - CrossSection xsec, CrossSection xsecErr, double evoptweight, double centralweight) const { - - //convert the event from the Herwig format to the HepMC format and write it to the common block - HepMC::GenEvent * ev = HepMCConverter::convert(*event, false,eUnit, lUnit); - - if(normoptweights) { evoptweight /= centralweight; } - - //reset the event - HepMCTraits::resetEvent(ev, no, evoptweight, event->optionalWeights()); - //set the cross section - HepMCTraits::setCrossSection(*ev,xsec/picobarn, - xsecErr/picobarn); - return ev; -} - -void FxFxAnalysis::analyze(ThePEG::tEventPtr event, long ieve, int loop, int state) { - Energy eUnit; - Length lUnit; - switch (_unitchoice) { - default: eUnit = GeV; lUnit = millimeter; break; - case 1: eUnit = MeV; lUnit = millimeter; break; - case 2: eUnit = GeV; lUnit = centimeter; break; - case 3: eUnit = MeV; lUnit = centimeter; break; - } - - tcFxFxEventHandlerPtr eh = dynamic_ptr_cast(event->primaryCollision()->handler()); - assert(eh); - - CrossSection xsec = eh->integratedXSec(); - CrossSection xsecErr = eh->integratedXSecErr(); - - optxsec = eh->optintegratedXSecMap(); - - - int ii = 0; - if(useoptweights) { - for (map::const_iterator it= optxsec.begin(); it!=optxsec.end(); ++it){ - OptXS[ii] = it->second/picobarn; - ii++; - } - } - - tSubProPtr sub = event->primarySubProcess(); - Ptr::tptr grp = - dynamic_ptr_cast::tptr>(sub); - - AnalysisHandler::analyze(event, ieve, loop, state); - // Rotate to CMS, extract final state particles and call analyze(particles). - // convert to hepmc - HepMC::GenEvent * hepmc = - makeEvent(event,sub,_nevent,eUnit,lUnit,xsec,xsecErr); - - //count weights here - std::vector strs; - if(_numweights == -999) { - _numweights = 0; - for (map::const_iterator it= event->optionalWeights().begin(); it!=event->optionalWeights().end(); ++it){ - string first_piece = it->first; - string word; - istringstream iss(first_piece, istringstream::in); - while( iss >> word ) strs.push_back(word); - if(strs[0] != "np") { _numweights++; } - strs.clear(); - } - } - double CentralWeight = 1.; - if(useoptweights) { - //find the weights - _i = 0;//counter - OptWeights.clear(); - for (map::const_iterator it= event->optionalWeights().begin(); it!=event->optionalWeights().end(); ++it){ - // std::cout << it->first << " => " << it->second << '\n'; - string first_piece = it->first; - string word; - istringstream iss(first_piece, istringstream::in); - while( iss >> word ) { - strs.push_back(word); - } - std::pair OptWeightsTemp; - /* if(strs[0] == "PDF") { - OptWeightsTemp.first = atoi(strs[2].c_str()); - OptWeightsTemp.second= it->second; - OptWeights.push_back(OptWeightsTemp); - _i++; - } - if(strs[0] == "SC") { - OptWeightsTemp.first = atoi(strs[3].c_str()); - OptWeightsTemp.second = it->second; - // cout << "OptWeightsTemp.first = " << OptWeightsTemp.first << " OptWeightsTemp.second = " << OptWeightsTemp.second << endl; - if(OptWeightsTemp.first == 1001) { cout << "OptWeightsTemp.second = " << OptWeightsTemp.second << endl; CentralWeight = OptWeightsTemp.second; } - OptWeights.push_back(OptWeightsTemp); - _i++; - }*/ - if(strs[0] != "np") { - OptWeightsTemp.first = atoi((it->first).c_str()); - OptWeightsTemp.second= it->second; - if(OptWeightsTemp.first == 1001) { /*cout << "OptWeightsTemp.second = " << OptWeightsTemp.second << endl;*/ CentralWeight = OptWeightsTemp.second; } - OptWeights.push_back(OptWeightsTemp); - _i++; - } - strs.clear(); - } - } - - /* multiple hepmcs for scale/pdf variations - */ - vector hepmcMULTI;// = new HepMC::GenEvent[_numweights]; - HepMC::GenEvent * hepmcMULTIi; - if(useoptweights) { - for(int rr = 0; rr < _numweights; rr++) { - double xsrr = optxsec[std::to_string(OptWeights[rr].first)]/picobarn; - /* cout << "xsec = " << xsec/picobarn << endl; - cout << "OptWeights[rr].second = " << OptWeights[rr].second << endl; - cout << "xsrr = " << xsrr << endl;*/ - hepmcMULTIi = makeEventW(event,_nevent,eUnit,lUnit,xsrr*picobarn,xsecErr,OptWeights[rr].second, CentralWeight); - hepmcMULTI.push_back(hepmcMULTIi); - } - } - - CurrentGenerator::Redirect stdout(cout); - - if ( _rivet ) _rivet->analyze(*hepmc); - - if(useoptweights) { - for(int rr = 0; rr < _numweights; rr++) { - if ( _rivetMULTI[rr] ) _rivetMULTI[rr]->analyze(*hepmcMULTI[rr]); - } - } - // delete hepmc events - delete hepmc; - if(useoptweights) { - for(int rr = 0; rr < _numweights; rr++) { - delete hepmcMULTI[rr]; - } - } - - if ( grp ) { - for ( SubProcessVector::const_iterator s = grp->dependent().begin(); - s != grp->dependent().end(); ++s ) { - - hepmc = makeEvent(event,*s,_nevent,eUnit,lUnit,xsec,xsecErr); - - if ( _rivet ) _rivet->analyze(*hepmc); - // delete hepmc event - delete hepmc; - - - } - - } - - ++_nevent; -} - - -ThePEG::IBPtr FxFxAnalysis::clone() const { - return new_ptr(*this); -} - -ThePEG::IBPtr FxFxAnalysis::fullclone() const { - return new_ptr(*this); -} - -void FxFxAnalysis::persistentOutput(ThePEG::PersistentOStream & os) const { - os << _analyses << filename << debug << useoptweights << normoptweights; -} - -void FxFxAnalysis::persistentInput(ThePEG::PersistentIStream & is, int) { - is >> _analyses >> filename >> debug >> useoptweights >> normoptweights; -} - -ThePEG::ClassDescription FxFxAnalysis::initFxFxAnalysis; -// Definition of the static class description member. - -void FxFxAnalysis::Init() { - static ThePEG::ClassDocumentation documentation - ("The FxFxAnalysis class is a simple class to allow analyses" - " from the Rivet library to be called from ThePEG"); - - static ThePEG::ParVector interfaceAnalyses - ("Analyses", - "The names of the Rivet analyses to use", - &FxFxAnalysis::_analyses, -1, "", "","" "", - false, false, ThePEG::Interface::nolimits); - - static Parameter interfaceRemnantId - ("RemnantId", - "Set the PDG id to be used for remnants.", - &FxFxAnalysis::_remnantId, 82, 0, 0, - false, false, Interface::nolimits); - - static Parameter interfaceFilename - ("Filename", - "The name of the file where the YODA histograms are put. If empty, " - "the run name will be used instead. '.yoda' will in any case be " - "appended to the file name.", - &FxFxAnalysis::filename, "", true, false); - - - static Switch interfaceDebug - ("Debug", - "Enable debug information from Rivet", - &FxFxAnalysis::debug, false, true, false); - static SwitchOption interfaceDebugNo - (interfaceDebug, - "No", - "Disable debug information.", - false); - static SwitchOption interfaceDebugYes - (interfaceDebug, - "Yes", - "Enable debug information from Rivet.", - true); - - - static Switch interfaceUseOptWeights - ("UseOptWeights", - "Enable debug information from Rivet", - &FxFxAnalysis::useoptweights, false, true, false); - static SwitchOption interfaceUseOptWeightsNo - (interfaceUseOptWeights, - "No", - "Disable optional weights", - false); - static SwitchOption interfaceUseOptWeightsYes - (interfaceUseOptWeights, - "Yes", - "Enable debug information from Rivet.", - true); - - static Switch interfaceNormOptWeights - ("NormOptWeights", - "Enable debug information from Rivet", - &FxFxAnalysis::normoptweights, false, true, false); - static SwitchOption interfaceNormOptWeightsNo - (interfaceNormOptWeights, - "No", - "Do not normalize optional weights", - false); - static SwitchOption interfacedNormOptWeightsYes - (interfaceNormOptWeights, - "Yes", - "Normalize optional weights.", - true); - - - - - interfaceAnalyses.rank(10); - -} - -void FxFxAnalysis::dofinish() { - AnalysisHandler::dofinish(); - if( _nevent > 0 && _rivet ) { - CurrentGenerator::Redirect stdout(cout); - _rivet->setCrossSection(generator()->integratedXSec()/picobarn); - _rivet->finalize(); - - string fname = filename; - fname = generator()->runName() + ".yoda"; - _rivet->writeData(fname); - } - - delete _rivet; - _rivet = 0; - - if(useoptweights) { - for(int rr = 0; rr < _numweights; rr++) { - cout << (OptWeights[rr].first) << ", cross section = " << OptXS[rr] << endl; - if( _nevent > 0 && _rivetMULTI[rr] ) { - double xsrr = optxsec[std::to_string(OptWeights[rr].first)]/picobarn; - // _rivetMULTI[rr]->setCrossSection(OptXS[rr]); - _rivetMULTI[rr]->setCrossSection(xsrr); - _rivetMULTI[rr]->finalize(); - - string fname = filename; - fname = generator()->runName() + "_" + std::to_string(OptWeights[rr].first) + ".yoda"; - - _rivetMULTI[rr]->writeData(fname); - } - delete _rivetMULTI[rr]; - _rivetMULTI[rr] = 0; - } - } -} - -void FxFxAnalysis::doinit() { - _numweights = -999; - - for(int rr = 0; rr < 120; rr++) { - OptXS.push_back(0.); - } - AnalysisHandler::doinit(); - if(_analyses.empty()) - throw ThePEG::Exception() << "Must have at least one analysis loaded in " - << "FxFxAnalysis::doinitrun()" - << ThePEG::Exception::runerror; - - // check that analysis list is available - _rivet = new Rivet::AnalysisHandler; //(fname); - _rivet->addAnalyses(_analyses); - if ( _rivet->analysisNames().size() != _analyses.size() ) { - throw ThePEG::Exception() - << "Rivet could not find all requested analyses.\n" - << "Use 'rivet --list-analyses' to check availability.\n" - << ThePEG::Exception::runerror; - } - - delete _rivet; - _rivet = 0; -} - -void FxFxAnalysis::doinitrun() { - _numweights = -999; - AnalysisHandler::doinitrun(); - // create FxFx analysis handler - if(useoptweights) { cout << "Warning: Using optional weights launches multiple rivet analyses. This may slow down your run substantially!" << endl; } - CurrentGenerator::Redirect stdout(cout); - _rivet = new Rivet::AnalysisHandler; //(fname); - _rivet->addAnalyses(_analyses); - if(useoptweights) { - for(int rr = 0; rr < 110; rr++) { - OptXS.push_back(0.); - _rivetMULTI[rr] = new Rivet::AnalysisHandler; //(fname); - _rivetMULTI[rr]->addAnalyses(_analyses); - } - } - - // check that analysis list is still available - if ( _rivet->analysisNames().size() != _analyses.size() ) { - throw ThePEG::Exception() - << "Rivet could not find all requested analyses.\n" - << "Use 'rivet --list-analyses' to check availability.\n" - << ThePEG::Exception::runerror; - } - if ( debug ) - Rivet::Log::setLevel("Rivet",Rivet::Log::DEBUG); -} diff --git a/MatrixElement/FxFx/FxFxAnalysis.h b/MatrixElement/FxFx/FxFxAnalysis.h deleted file mode 100644 --- a/MatrixElement/FxFx/FxFxAnalysis.h +++ /dev/null @@ -1,262 +0,0 @@ -// -*- C++ -*- -#ifndef THEPEG_FxFxAnalysis_H -#define THEPEG_FxFxAnalysis_H -// -// This is the declaration of the FxFxAnalysis class. -// - -#include "ThePEG/Handlers/AnalysisHandler.h" -#include "Rivet/AnalysisHandler.hh" -#include "FxFxReader.h" -#include "FxFxEventHandler.h" - -namespace ThePEG { - -/** - * Here is the documentation of the FxFxAnalysis class. - * - * @see \ref FxFxAnalysisInterfaces "The interfaces" - * defined for FxFxAnalysis. - */ -class FxFxAnalysis: public ThePEG::AnalysisHandler { - -public: - - /** - * The default constructor. - */ - FxFxAnalysis(); - -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(ThePEG::tEventPtr event, long ieve, int loop, int state); - - /** - * Produca a HepMC event for the given subprocess - */ - HepMC::GenEvent * makeEvent(tEventPtr event, tSubProPtr sub, long no, - Energy eUnit, Length lUnit, - CrossSection xsec, CrossSection xsecErr) const; - - /** - * Produca a HepMC event for the given subprocess - */ - HepMC::GenEvent * makeEventW(tEventPtr event, long no, - Energy eUnit, Length lUnit, - CrossSection xsec, CrossSection xsecErr, double evoptweight, double centralweight) const; - - int _i; - - int _numweights; - - //@} - -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(ThePEG::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(ThePEG::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. - */ - virtual ThePEG::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. - */ - virtual ThePEG::IBPtr fullclone() const; - //@} - -protected: - - /** @name Standard Interfaced functions. */ - //@{ - /** - * Initialize this object. Called in the read phase. - */ - 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(); - //@} - -private: - - /** - * The static object used to initialize the description of this class. - * Indicates that this is a concrete class with persistent data. - */ - static ThePEG::ClassDescription initFxFxAnalysis; - - /** - * The assignment operator is private and must never be called. - * In fact, it should not even be implemented. - */ - FxFxAnalysis & operator=(const FxFxAnalysis &) = delete; - -private: - /** - * The PDG ID to be used for remnants - */ - long _remnantId; - - /** - * The HepMC format - */ - int _format; - - /** - * Selector for the choice of units - */ - int _unitchoice; - - /** - * Choice of output precision in GenEvent format - */ - unsigned int _geneventPrecision; - - /** - * The Analyses to use - */ - vector _analyses; - - /** - * The base name of the output file. - */ - string filename; - - /** - * Enable debugging information from FxFx - */ - bool debug; - - /** - * Enable use of optional weights in analysis - */ - bool useoptweights = false; - - /** - * normalize optional weights to the central weight - */ - - bool normoptweights = false; - - /** - * The FxFxAnalysisHandler - */ - Rivet::AnalysisHandler * _rivet; - - /** - * The FxFxAnalysisHandlers for multiple weights - */ - Rivet::AnalysisHandler * _rivetMULTI[120]; - - - /** - * holders of weights and cross section - */ - - std::vector< std::pair > OptWeights; - std::vector OptXS; - map optxsec; - - - /** - * Event count - */ - unsigned long _nevent; -}; - -} - -#include "ThePEG/Utilities/ClassTraits.h" - -namespace ThePEG { - -/** @cond TRAITSPECIALIZATIONS */ - -/** This template specialization informs ThePEG about the - * base classes of FxFxAnalysis. */ -template <> -struct BaseClassTrait { - /** Typedef of the first base class of FxFxAnalysis. */ - typedef AnalysisHandler NthBase; -}; - -/** This template specialization informs ThePEG about the name of - * the FxFxAnalysis 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::FxFxAnalysis"; } - /** - * The name of a file containing the dynamic library where the class - * FxFxAnalysis is implemented. It may also include several, space-separated, - * libraries if the class FxFxAnalysis depends on other classes (base classes - * excepted). In this case the listed libraries will be dynamically - * linked in the order they are specified. - */ - static string library() { return "FxFxAnalysis.so"; } -}; - -/** @endcond */ - -} - -#endif /* THEPEG_FxFxAnalysis_H */