#include "SHERPA/Main/Sherpa.H" #include "SHERPA/Initialization/Initialization_Handler.H" #include "SHERPA/Single_Events/Event_Handler.H" #include "SHERPA/Single_Events/Analysis_Phase.H" #include "SHERPA/Single_Events/Userhook_Phase.H" #include "SHERPA/Single_Events/Output_Phase.H" #include "SHERPA/Single_Events/EvtReadin_Phase.H" #include "SHERPA/Single_Events/Signal_Processes.H" #include "SHERPA/Single_Events/Hard_Decays.H" #include "SHERPA/Single_Events/Minimum_Bias.H" #include "SHERPA/Single_Events/Multiple_Interactions.H" #include "SHERPA/Single_Events/Jet_Evolution.H" #include "SHERPA/Single_Events/Signal_Process_FS_QED_Correction.H" #include "SHERPA/Single_Events/Beam_Remnants.H" #include "SHERPA/Single_Events/Hadronization.H" #include "SHERPA/Single_Events/Hadron_Decays.H" #include "SHERPA/PerturbativePhysics/Hard_Decay_Handler.H" #include "SHERPA/Tools/HepMC2_Interface.H" #include "SHERPA/Tools/HepMC3_Interface.H" #include "PHASIC++/Decays/Decay_Channel.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Shell_Tools.H" #include "ATOOLS/Org/Library_Loader.H" #include "ATOOLS/Org/My_MPI.H" #include "ATOOLS/Org/CXXFLAGS.H" #include "ATOOLS/Org/CXXFLAGS_PACKAGES.H" #include "ATOOLS/Org/Scoped_Settings.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Phys/KF_Table.H" #include "PDF/Main/ISR_Handler.H" #include "PDF/Main/PDF_Base.H" #include <cstring> using namespace SHERPA; using namespace ATOOLS; using namespace std; Sherpa::Sherpa(int argc, char* argv[]) : p_inithandler(nullptr), p_eventhandler(nullptr) #ifdef USING__HEPMC2 , p_hepmc2(nullptr) #endif #ifdef USING__HEPMC3 , p_hepmc3(nullptr) #endif { ATOOLS::mpi = new My_MPI(); ATOOLS::exh = new Terminator_Object_Handler(); ATOOLS::msg = new Message(); // rpa should be constructed before initializing the main settings, since the // latter might throw an exception and rpa would be involved in terminating // the program then; however, do not call its Init method yet, because this // in turn needs the Settings to be initialized ATOOLS::rpa = new Run_Parameter(); Settings::InitializeMainSettings(argc, argv); ATOOLS::ran = new Random(1234); ATOOLS::s_loader = new Library_Loader(); PDF::pdfdefs = new PDF::PDF_Defaults(); m_trials = 0; m_debuginterval = 0; m_debugstep = -1; m_displayinterval = 100; m_evt_starttime = -1.0; exh->AddTerminatorObject(this); } Sherpa::~Sherpa() { if (msg_LevelIsInfo()) { Return_Value::PrintStatistics(msg->Out()); if (p_inithandler->GetVariations()) { p_inithandler->GetVariations()->PrintStatistics(msg->Out()); } } Blob_List::PrintMomFailStatistics(msg->Out()); msg->PrintRates(); PHASIC::Decay_Channel::PrintMaxKinFailStatistics(msg->Out()); rpa->gen.WriteCitationInfo(); if (p_eventhandler) { delete p_eventhandler; p_eventhandler = nullptr; } if (p_inithandler) { delete p_inithandler; p_inithandler = nullptr; } #ifdef USING__HEPMC2 if (p_hepmc2) { delete p_hepmc2; p_hepmc2 = nullptr; } #endif #ifdef USING__HEPMC3 if (p_hepmc3) { delete p_hepmc3; p_hepmc3 = NULL; } #endif Settings& s = Settings::GetMainSettings(); if (s["CHECK_SETTINGS"].SetDefault(true).Get<bool>()) Settings::FinalizeMainSettings(); exh->RemoveTerminatorObject(this); delete ATOOLS::s_loader; delete PDF::pdfdefs; delete ATOOLS::rpa; delete ATOOLS::ran; #ifdef USING__MPI mpi->Barrier(); #endif delete ATOOLS::msg; delete ATOOLS::exh; delete ATOOLS::mpi; for (KF_Table::const_iterator kfit(s_kftable.begin());kfit!=s_kftable.end();++kfit) delete kfit->second; ATOOLS::s_kftable.clear(); } bool Sherpa::InitializeTheRun() { Settings& s = Settings::GetMainSettings(); p_inithandler = new Initialization_Handler(); RegisterDefaults(); mpi->PrintRankInfo(); DrawLogo(s["PRINT_VERSION_INFO"].Get<bool>()); int initonly=s["INIT_ONLY"].Get<int>(); if (initonly) rpa->gen.SetNumberOfEvents(0); if (p_inithandler->InitializeTheFramework()) { if (initonly==1) THROW(normal_exit,"Initialization complete."); if (initonly==2) return true; if (!p_inithandler->CalculateTheHardProcesses()) return false; m_showtrials=s["SHOW_NTRIALS"].Get<bool>(); // read in from status path bool res(true); std::string statuspath(s["STATUS_PATH"].Get<std::string>()); if (statuspath != "") { res=exh->ReadInStatus(statuspath); } m_debuginterval = s["DEBUG_INTERVAL"].Get<long int>(); m_debugstep = s["DEBUG_STEP"].Get<long int>(); m_displayinterval=s["EVENT_DISPLAY_INTERVAL"].Get<int>(); m_evt_output = s["EVT_OUTPUT"].Get<int>(); m_evt_output_start = s["EVT_OUTPUT_START"].Get<int>(); return res; } msg_Error()<<"Error in Sherpa::InitializeRun()"<<endl <<" Did not manage to initialize the framework."<<endl <<" Try to run nevertheless ... ."<<endl; return 0; } void Sherpa::RegisterDefaults() { Settings& s = Settings::GetMainSettings(); s["PRINT_VERSION_INFO"].SetDefault(false); s["INIT_ONLY"].SetDefault(0); s["SHOW_NTRIALS"].SetDefault(false); s["DEBUG_INTERVAL"].SetDefault(0); s["DEBUG_STEP"].SetDefault(-1); s["EVENT_DISPLAY_INTERVAL"].SetDefault(100); s["EVT_OUTPUT"].SetDefault(msg->Level()); s["MSG_LIMIT"].SetDefault(20); msg->SetLimit(s["MSG_LIMIT"].Get<int>()); const int evtoutput{ s["EVT_OUTPUT"].Get<int>() }; s["EVT_OUTPUT_START"].SetDefault(evtoutput != msg->Level() ? 1 : 0); } bool Sherpa::InitializeTheEventHandler() { eventtype::code mode = p_inithandler->Mode(); p_eventhandler = new Event_Handler(); Analysis_Vector *anas(p_inithandler->GetAnalyses()); for (Analysis_Vector::iterator it=anas->begin(); it!=anas->end(); ++it) { (*it)->SetEventHandler(p_eventhandler); } if (mode==eventtype::EventReader) { p_eventhandler->AddEventPhase(new EvtReadin_Phase(p_inithandler->GetEventReader())); p_eventhandler->AddEventPhase(new Hard_Decays(p_inithandler->GetHardDecayHandler())); p_eventhandler->AddEventPhase(new Beam_Remnants(p_inithandler->GetBeamRemnantHandler())); } else { p_eventhandler->AddEventPhase(new Signal_Processes(p_inithandler->GetMatrixElementHandler())); p_eventhandler->AddEventPhase(new Minimum_Bias(p_inithandler->GetSoftCollisionHandlers())); p_eventhandler->AddEventPhase(new Hard_Decays(p_inithandler->GetHardDecayHandler())); p_eventhandler->AddEventPhase(new Jet_Evolution(p_inithandler->GetMatrixElementHandler(), p_inithandler->GetHardDecayHandler(), p_inithandler->GetHDHandler(), p_inithandler->GetMIHandlers(), p_inithandler->GetSoftCollisionHandlers(), p_inithandler->GetShowerHandlers(), p_inithandler->GetRemnantHandlers())); p_eventhandler->AddEventPhase(new Signal_Process_FS_QED_Correction( p_inithandler->GetMatrixElementHandler(), p_inithandler->GetSoftPhotonHandler())); p_eventhandler->AddEventPhase(new Multiple_Interactions(p_inithandler->GetMIHandlers())); p_eventhandler->AddEventPhase(new Beam_Remnants(p_inithandler->GetBeamRemnantHandler())); p_eventhandler->AddEventPhase(new Hadronization(p_inithandler->GetColourReconnectionHandler(), p_inithandler->GetFragmentation())); p_eventhandler->AddEventPhase(new Hadron_Decays(p_inithandler->GetHDHandler())); } p_eventhandler->AddEventPhase(new Userhook_Phase(this)); if (!anas->empty()) p_eventhandler->AddEventPhase(new Analysis_Phase(anas)); if (!p_inithandler->GetOutputs()->empty()) p_eventhandler->AddEventPhase(new Output_Phase(p_inithandler->GetOutputs(), p_eventhandler)); p_eventhandler->SetFilter(p_inithandler->GetFilter()); p_eventhandler->PrintGenericEventStructure(); ran->EraseLastIncrementedSeed(); return 1; } bool Sherpa::GenerateOneEvent(bool reset) { if (m_evt_output_start>0 && m_evt_output_start==rpa->gen.NumberOfGeneratedEvents()+1) { msg->SetLevel(m_evt_output); } if(m_debuginterval>0 && rpa->gen.NumberOfGeneratedEvents()%m_debuginterval==0 && (p_inithandler->GetMatrixElementHandler()->SeedMode()!=3 || rpa->gen.NumberOfGeneratedEvents()==0)) { std::string fname=ToString(rpa->gen.NumberOfGeneratedEvents())+".dat"; ran->WriteOutStatus(("random."+fname).c_str()); } if (m_debugstep>=0) { if (p_inithandler->GetMatrixElementHandler()->SeedMode()!=3) ran->ReadInStatus(("random."+ToString(m_debugstep)+".dat").c_str()); else { ran->ReadInStatus("random.0.dat"); ran->FastForward(m_debugstep); } } if (m_evt_starttime<0.0) m_evt_starttime=rpa->gen.Timer().RealTime(); if (reset) p_eventhandler->Reset(); if (p_eventhandler->GenerateEvent(p_inithandler->Mode())) { if(m_debuginterval>0 && rpa->gen.NumberOfGeneratedEvents()%m_debuginterval==0){ std::string fname=ToString(rpa->gen.NumberOfGeneratedEvents())+".dat"; std::ofstream eventout(("refevent."+fname).c_str()); eventout<<"# trial "<<rpa->gen.NumberOfTrials()-1<<std::endl; eventout<<*p_eventhandler->GetBlobs()<<std::endl; eventout.close(); } if (m_debugstep>=0) { std::ofstream event(("event."+ToString(m_debugstep)+".dat").c_str()); event<<*p_eventhandler->GetBlobs()<<std::endl; event.close(); THROW(normal_exit,"Debug event written."); } rpa->gen.SetNumberOfGeneratedEvents(rpa->gen.NumberOfGeneratedEvents()+1); Blob_List *blobs(p_eventhandler->GetBlobs()); /// Increase m_trials --- based on signal blob["Trials"] if existent if (blobs->FindFirst(btp::Signal_Process) == nullptr) { m_trials+=1; msg_Debugging()<<" No Signal_Process Blob found, increasing m_trials by 1\n"; } else { m_trials+=(*blobs->FindFirst(btp::Signal_Process))["Trials"]->Get<double>(); } if (msg_LevelIsEvents()) { if (!blobs->empty()) { msg_Out()<<" -------------------------------------------------\n"; for (Blob_List::iterator blit=blobs->begin(); blit!=blobs->end();++blit) msg_Out()<<*(*blit)<<std::endl; msg_Out()<<" -------------------------------------------------\n"; } else msg_Out()<<" ******** Empty event ******** "<<std::endl; } int i=rpa->gen.NumberOfGeneratedEvents(); int nevt=rpa->gen.NumberOfEvents(); msg_Events()<<"Sherpa : Passed "<<i<<" events."<<std::endl; int exp; for (exp=5; i/int(pow(10,exp))==0; --exp) {} if (((rpa->gen.BatchMode()&4 && i%m_displayinterval==0) || (!(rpa->gen.BatchMode()&4) && i%int(pow(10,exp))==0)) && i<rpa->gen.NumberOfEvents()) { double diff=rpa->gen.Timer().RealTime()-m_evt_starttime; msg_Info()<<" Event "<<i; if (m_showtrials) msg_Info()<<"("+ToString(m_trials)+")"; msg_Info()<<" ( "; if (rpa->gen.BatchMode()&16) { msg_Info()<<diff<<"s elapsed / " <<((nevt-i)/(double)i*diff)<<"s"; } else { msg_Info()<<FormatTime(size_t(diff))<<" elapsed / " <<FormatTime(size_t((nevt-i)/(double)i*diff)); } msg_Info()<<" left ) -> ETA: "<<rpa->gen.Timer(). StrFTime("%a %b %d %H:%M",time_t((nevt-i)/(double)i*diff))<<" "; p_eventhandler->PerformMemoryMonitoring(); const Uncertain<double> xs = p_eventhandler->TotalNominalXSMPI(); if (!(rpa->gen.BatchMode()&2)) msg_Info()<<"\n "; msg_Info() << "XS = " << xs.value << " pb +- ( " << xs.error << " pb = " << xs.PercentError() << " % ) "; if (rpa->gen.BatchMode()&8) msg_Info()<<" Process was "<<p_eventhandler->CurrentProcess()<<" "; if (!(rpa->gen.BatchMode()&2)) msg_Info()<<mm(1,mm::up); if (rpa->gen.BatchMode()&2) { msg_Info()<<std::endl; } else { msg_Info()<<bm::cr<<std::flush; } } return 1; } return 0; } #ifdef USING__HEPMC2 void Sherpa::FillHepMCEvent(HepMC::GenEvent& event) { if (!p_hepmc2) p_hepmc2 = new SHERPA::HepMC2_Interface(); ATOOLS::Blob_List* blobs = GetEventHandler()->GetBlobs(); p_hepmc2->Sherpa2HepMC(blobs, event); p_hepmc2->AddCrossSection(event, p_eventhandler->TotalNominalXS()); } #endif #ifdef USING__HEPMC3 void Sherpa::FillHepMCEvent(HepMC3::GenEvent& event) { if (p_hepmc3==NULL) p_hepmc3 = new SHERPA::HepMC3_Interface(); ATOOLS::Blob_List* blobs=GetEventHandler()->GetBlobs(); p_hepmc3->Sherpa2HepMC(blobs, event); p_hepmc3->AddCrossSection(event, p_eventhandler->TotalXS(), p_eventhandler->TotalErr()); } #endif Uncertain<double> Sherpa::TotalNominalXS() const { return p_eventhandler->TotalNominalXS(); } std::string Sherpa::PDFInfo() { std::string pdf="Unknown"; PDF::ISR_Handler* isr=GetInitHandler()->GetISRHandler(PDF::isr::hard_process); if (isr) { if (isr->PDF(0)) { pdf=isr->PDF(0)->Type(); if (isr->PDF(1) && isr->PDF(1)->Type()!=pdf) { pdf="Unknown"; } } } return pdf; } void Sherpa::PrepareTerminate() { SummarizeRun(); exh->RemoveTerminatorObject(this); } bool Sherpa::SummarizeRun() { if (p_eventhandler) { msg_Info()<<" Event "<<rpa->gen.NumberOfGeneratedEvents()<<" ( " <<size_t(rpa->gen.Timer().RealTime()-m_evt_starttime) <<" s total ) = " << rpa->gen.NumberOfGeneratedEvents()*3600*24/ ((size_t) rpa->gen.Timer().RealTime()-m_evt_starttime) <<" evts/day "<<std::endl; p_eventhandler->Finish(); } return true; } long int Sherpa::NumberOfEvents() const { return rpa->gen.NumberOfEvents(); } const Blob_List &Sherpa::GetBlobList() const { return *p_eventhandler->GetBlobs(); } double Sherpa::GetMEWeight(const Cluster_Amplitude &l,const int mode) const { return p_inithandler->GetMatrixElementHandler()-> GetWeight(ampl,ATOOLS::nlo_type::lo,mode); } void Sherpa::DrawLogo(const bool& shouldprintversioninfo) { msg_Info()<<"-----------------------------------------------------------------------------"<<std::endl; if (msg->Level()>0) msg_Out()<<"----------- Event generation run with SHERPA started ....... -----------"<<std::endl; msg_Info()<<"-----------------------------------------------------------------------------"<<std::endl <<"................................................ | + "<<std::endl <<"................................................ || | + + "<<std::endl <<"................................... .... | | / + "<<std::endl <<"................. ................ _,_ | .... || +| + + "<<std::endl <<"............................... __.' ,\\| ... || / +| + "<<std::endl <<".............................. ( \\ \\ ... | | | + + \\ + "<<std::endl <<"............................. ( \\ -/ .... || + | + "<<std::endl <<"........ ................... <S /()))))~~~~~~~~## + /\\ + "<<std::endl <<"............................ (!H (~~)))))~~~~~~#/ + + | + "<<std::endl <<"................ ........... (!E (~~~))))) /|/ + + "<<std::endl <<"............................ (!R (~~~))))) ||| + + + "<<std::endl <<"..... ...................... (!P (~~~~))) /| + + + "<<std::endl <<"............................ (!A> (~~~~~~~~~## + + + "<<std::endl <<"............................. ~~(! '~~~~~~~ \\ + + + + "<<std::endl <<"............................... `~~~QQQQQDb // | + + + + "<<std::endl <<"........................ .......... IDDDDP|| \\ + + + + + +"<<std::endl <<".................................... IDDDI|| \\ + "<<std::endl <<".................................... IHD HD|| \\ + + + + + + + +"<<std::endl <<"................................... IHD ##| :-) + +\\ + "<<std::endl <<"......... ............... ......... IHI ## / / + + + + +\\ + "<<std::endl <<"................................... IHI/ / / + + + + + "<<std::endl <<"................................... ## | | / / + + + + / + "<<std::endl <<"....................... /TT\\ ..... ##/ /// / + + + + + + +/ + "<<std::endl <<"......................./TTT/T\\ ... /TT\\/\\\\\\ / + + + + + + +/ \\ + "<<std::endl <<"....................../TTT/TTTT\\...|TT/T\\\\\\/ + ++ + / "<<std::endl <<"-----------------------------------------------------------------------------"<<std::endl <<std::endl <<" SHERPA version "<<SHERPA_VERSION<<"."<<SHERPA_SUBVERSION<<" ("<<SHERPA_NAME<<")"<<std::endl <<" "<<std::endl <<" Authors: Enrico Bothmann, Stefan Hoeche, Frank Krauss, "<<std::endl <<" Silvan Kuttimalai, Marek Schoenherr, Holger Schulz, "<<std::endl <<" Steffen Schumann, Frank Siegert, Korinna Zapp "<<std::endl <<" Former Authors: Timo Fischer, Tanju Gleisberg, Hendrik Hoeth, "<<std::endl <<" Ralf Kuhn, Thomas Laubrich, Andreas Schaelicke, "<<std::endl <<" Jan Winter "<<std::endl <<" "<<std::endl <<" This program uses a lot of genuine and original research work "<<std::endl <<" by other people. Users are encouraged to refer to "<<std::endl <<" the various original publications. "<<std::endl <<" "<<std::endl <<" Users are kindly asked to refer to the documentation "<<std::endl <<" published under JHEP 02(2009)007 "<<std::endl <<" "<<std::endl <<" Please visit also our homepage "<<std::endl <<" "<<std::endl <<" http://sherpa.hepforge.org "<<std::endl <<" "<<std::endl <<" for news, bugreports, updates and new releases. "<<std::endl <<" "<<std::endl <<"-----------------------------------------------------------------------------"<<std::endl; rpa->gen.PrintGitVersion(msg->Info(), shouldprintversioninfo); rpa->gen.AddCitation (0,"The complete Sherpa package is published under \\cite{Gleisberg:2008ta}."); }