Commit 0eebe61f authored by Stefan Hoeche's avatar Stefan Hoeche

svn merge -r30895:31449 ^/branches/ntuple-variations .

parent 23b23ec1
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Decays","trunk/SHERPA","31416M","a996e1db8b65efc861a3d529d779a078");
("AHADIC++/Decays","trunk/SHERPA","31449M","a996e1db8b65efc861a3d529d779a078");
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Formation","trunk/SHERPA","31416M","ab96bb66a9ab2455d9e33233eeee4ab7");
("AHADIC++/Formation","trunk/SHERPA","31449M","4151c714418bf46c6e048a04fe96f73c");
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Main","trunk/SHERPA","31416M","6aec0d693c797d36f254e1137a42f6f9");
("AHADIC++/Main","trunk/SHERPA","31449M","6aec0d693c797d36f254e1137a42f6f9");
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Tools","trunk/SHERPA","31416M","96fdcbe6b883bea5fee2f27a7ff54f2e");
("AHADIC++/Tools","trunk/SHERPA","31449M","017ab3ca355adea8fcf06f3ca8f99506");
......@@ -170,6 +170,9 @@ bool Analysis_Handler::Init()
mode=mode|m_weighted;
m_analyses.push_back(new Primitive_Analysis(this,ToString(i),mode));
m_analyses.back()->SetOutputPath(outpath);
int usedb;
if (!reader.ReadFromFile(usedb,"USE_DB")) usedb=0;
m_analyses.back()->SetUseDB(usedb);
std::string maxjettag;
if (!reader.ReadFromFile(maxjettag,"NMAX_JETS")) maxjettag="";
m_analyses.back()->SetMaxJetTag(maxjettag);
......@@ -263,7 +266,7 @@ bool Analysis_Handler::WriteOut()
}
for (Analyses_Vector::const_iterator ait=m_analyses.begin();
ait!=m_analyses.end();++ait) {
(*ait)->FinishAnalysis(OutputPath());
(*ait)->FinishAnalysis(OutputPath(),0);
(*ait)->RestoreAnalysis();
}
return true;
......@@ -286,7 +289,7 @@ bool Analysis_Handler::Finish()
ait!=m_analyses.end();++ait) {
msg_Info()<<" Writing to '"<<OutputPath()<<(*ait)->OutputPath()
<<"'."<<std::endl;
(*ait)->FinishAnalysis(OutputPath());
(*ait)->FinishAnalysis(OutputPath(),0);
(*ait)->RestoreAnalysis();
}
msg_Info()<<"}"<<std::endl;
......
......@@ -25,6 +25,7 @@ Primitive_Analysis::Primitive_Analysis
m_nevt = 0;
p_partner = this;
m_mode = mode;
m_usedb = 0;
m_name = std::string("Analysis : ") + _name;
msg_Tracking()<<" Initializing Primitive_Analysis : "<<m_name<<std::endl;
......@@ -35,6 +36,7 @@ Primitive_Analysis::Primitive_Analysis(Analysis_Handler *const ana,const int mod
{
p_ana=ana;
m_mode = mode;
m_usedb = 0;
m_name = std::string("Analysis : noname");
msg_Tracking()<<" Initializing Primitive_Analysis : "<<m_name<<std::endl;
......@@ -91,6 +93,7 @@ Primitive_Analysis * Primitive_Analysis::GetSubAnalysis
Primitive_Analysis * ana = new Primitive_Analysis(p_ana,m_name.substr(11)+key,mode);
if (master) ana->SetPartner(p_partner);
ana->SetMaxJetTag(m_maxjettag);
ana->SetVarId(m_varid);
for (size_t i=0;i<m_objects.size();i++) {
if (m_objects[i]->IsObservable() || !master)
......@@ -172,8 +175,6 @@ void Primitive_Analysis::CallSubAnalysis(const Blob_List * const bl, double valu
mode=m_mode^ANALYSIS::splitt_jetseeds;
if (!m_splitjetconts)
mode=mode-(mode&ANALYSIS::output_this);
if (name.find("__")==std::string::npos) key="X";
else {
std::string fsname(name.substr(name.find("__")+3));
fsname=fsname.substr(fsname.find("__")+3);
fsname=fsname.substr(fsname.find("__")+2);
......@@ -184,7 +185,6 @@ void Primitive_Analysis::CallSubAnalysis(const Blob_List * const bl, double valu
key=JetID(fsname,m_maxjettag);
key="j"+key;
}
}
else {
mode=m_mode^ANALYSIS::splitt_process;
// if (m_mode&ANALYSIS::output_process) mode=mode|ANALYSIS::output_this;
......@@ -202,22 +202,21 @@ void Primitive_Analysis::CallSubAnalysis(const Blob_List * const bl, double valu
}
void Primitive_Analysis::DoAnalysis(const Blob_List * const bl, double weight)
void Primitive_Analysis::DoAnalysis(const Blob_List * const bl, const double value)
{
p_sub=p_real=NULL;
++m_nevt;
m_called.clear();
if (IsNan(weight)) {
if (IsNan(value)) {
msg_Error()<<METHOD<<"(): Event weight is nan. Skip."<<std::endl;
return;
}
if (m_mode&ANALYSIS::split_vars) {
m_mode&=~ANALYSIS::split_vars;
Blob *sp(bl->FindFirst(btp::Signal_Process));
Blob_Data_Base *info((*sp)["Variation_Weights"]);
if (info) {
int mode=m_mode|ANALYSIS::output_this;
int mode=(m_mode^ANALYSIS::split_vars)|ANALYSIS::output_this;
ATOOLS::Variation_Weights vars(info->Get<Variation_Weights>());
m_nvar=vars.GetNumberOfVariations();
if (m_nvar) {
......@@ -225,7 +224,10 @@ void Primitive_Analysis::DoAnalysis(const Blob_List * const bl, double weight)
std::string name(vars.GetVariationNameAt(i));
Primitive_Analysis *ana=GetSubAnalysis(bl,name,mode,false);
ana->SetVarId(i);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
return;
}
}
}
......@@ -234,41 +236,33 @@ void Primitive_Analysis::DoAnalysis(const Blob_List * const bl, double weight)
int mode=m_mode^ANALYSIS::splitt_phase;
if (m_mode&ANALYSIS::do_me) {
Primitive_Analysis *ana(GetSubAnalysis(bl,"ME",mode));
ana->DoAnalysis(bl,weight);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
if (m_mode&ANALYSIS::do_menlo) {
Primitive_Analysis *ana(GetSubAnalysis(bl,"MENLO",mode));
ana->DoAnalysis(bl,weight);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
if (m_mode&ANALYSIS::do_mi) {
Primitive_Analysis *ana(GetSubAnalysis(bl,"MI",mode));
ana->DoAnalysis(bl,weight);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
if (m_mode&ANALYSIS::do_shower) {
Primitive_Analysis *ana(GetSubAnalysis(bl,"Shower",mode));
ana->DoAnalysis(bl,weight);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
if (m_mode&ANALYSIS::do_hadron) {
Primitive_Analysis *ana(GetSubAnalysis(bl,"Hadron",mode));
ana->DoAnalysis(bl,weight);
ana->DoAnalysis(bl,value);
m_called.insert(ana);
}
return;
}
if (m_mode&ANALYSIS::do_menlo) {
if (DoAnalysisNLO(bl,weight)) {
for (Analysis_List::iterator it=m_subanalyses.begin();
it!=m_subanalyses.end();++it)
if (m_varid==-1 && it->second->VarId()>-1) {
it->second->DoAnalysis(bl,weight);
m_called.insert(it->second);
}
return;
}
if (DoAnalysisNLO(bl,value)) return;
}
ClearAllData();
......@@ -287,8 +281,9 @@ void Primitive_Analysis::DoAnalysis(const Blob_List * const bl, double weight)
Blob *sp(bl->FindFirst(btp::Signal_Process));
// if no signal process present (i.e. hadrons execs etc.),
// assume weight=1, ncount=1
double ncount(1.);
double weight(1.), ncount(1.);
if (sp) {
weight=(*sp)["Weight"]->Get<double>();
ncount=(*sp)["Trials"]->Get<double>();
}
if (m_varid>-1) {
......@@ -410,8 +405,9 @@ bool Primitive_Analysis::DoAnalysisNLO(const Blob_List * const bl, const double
}
void Primitive_Analysis::FinishAnalysis(const std::string & resdir)
void Primitive_Analysis::FinishAnalysis(const std::string & resdir,int mode)
{
if (m_usedb) mode=1;
#ifdef USING__MPI
if (MPI::COMM_WORLD.Get_size()>1) {
int rank=MPI::COMM_WORLD.Get_rank();
......@@ -440,34 +436,55 @@ void Primitive_Analysis::FinishAnalysis(const std::string & resdir)
}
if (MPI::COMM_WORLD.Get_rank()==0)
#endif
ATOOLS::MakeDir(resdir+OutputPath());
{
if (m_usedb) {
My_In_File::OpenDB(resdir+OutputPath()+"/");
My_In_File::ExecDB(resdir+OutputPath(),"begin");
}
if (mode==0) ATOOLS::MakeDir(resdir+OutputPath());
}
if (m_mode&ANALYSIS::do_menlo) {
if ((m_mode&ANALYSIS::split_vars) && m_nvar) {
for (Analysis_List::iterator it=m_subanalyses.begin();
it!=m_subanalyses.end();++it)
if (m_varid==-1 && it->second->VarId()>-1) {
it!=m_subanalyses.end();++it) {
std::string dir=resdir+OutputPath()+std::string("/")+it->first;
it->second->FinishAnalysis(dir);
it->second->FinishAnalysis(dir,mode);
}
if (m_usedb) {
My_In_File::ExecDB(resdir+OutputPath(),"commit");
My_In_File::CloseDB(resdir+OutputPath()+"/");
}
return;
}
for (size_t i=0;i<m_objects.size();i++) {
m_objects[i]->EndEvaluation(1.);
m_objects[i]->Output(resdir+OutputPath());
}
if (m_usedb) {
My_In_File::ExecDB(resdir+OutputPath(),"commit");
My_In_File::CloseDB(resdir+OutputPath()+"/");
}
return;
}
for (Analysis_List::iterator it=m_subanalyses.begin();
it!=m_subanalyses.end();++it) {
std::string dir=resdir+OutputPath()+std::string("/")+it->first;
it->second->FinishAnalysis(dir);
it->second->FinishAnalysis(dir,mode);
}
if (!(m_mode&ANALYSIS::splitt_phase)) {
if (!(m_mode&ANALYSIS::splitt_phase) &&
!((m_mode&ANALYSIS::split_vars) && m_nvar)) {
for (size_t i=0;i<m_objects.size();i++) {
m_objects[i]->EndEvaluation();
if (m_mode&ANALYSIS::output_this)
m_objects[i]->Output(resdir+OutputPath());
}
}
if (m_usedb) {
My_In_File::ExecDB(resdir+OutputPath(),"commit");
My_In_File::CloseDB(resdir+OutputPath()+"/");
}
}
void Primitive_Analysis::RestoreAnalysis()
......
......@@ -51,7 +51,7 @@ namespace ANALYSIS {
class Primitive_Analysis: public ATOOLS::File_IO_Base {
private:
int m_mode, m_varid, m_nvar;
int m_mode, m_varid, m_nvar, m_usedb;
long m_nevt;
std::string m_name, m_maxjettag;
......@@ -95,7 +95,7 @@ namespace ANALYSIS {
void DoAnalysis(const ATOOLS::Blob_List * const , double);
bool DoAnalysisNLO(const ATOOLS::Blob_List * const , const double);
void FinishAnalysis(const std::string &);
void FinishAnalysis(const std::string &,int mode);
void RestoreAnalysis();
void AddObject(Analysis_Object *obj);
......@@ -135,6 +135,8 @@ namespace ANALYSIS {
inline const int &VarId() const { return m_varid; }
inline void SetUseDB(const int &usedb) { m_usedb=usedb; }
};// end of class Primitive_Analysis
}// end of namespace ANALYSIS
......
......@@ -5,6 +5,7 @@
#include <termios.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <vector>
using namespace ATOOLS;
......@@ -35,6 +36,16 @@ void PrintInfo()
#include "ATOOLS/Org/Exception.H"
#include "ATOOLS/Org/My_MPI.H"
#include "ATOOLS/Org/My_File.H"
#include <sqlite3.h>
int AddFiles(void *data,int argc,char **argv,char **name)
{
if (argc!=1 || strcmp(name[0],"file")) return 1;
msg_IODebugging()<<" '"<<argv[0]<<"'\n";
static_cast<std::vector<std::string>*>(data)->push_back(argv[0]);
return 0;
}
int main(int argc,char **argv)
{
......@@ -89,13 +100,45 @@ int main(int argc,char **argv)
i++;
vector<string> inlist;
for (;i<argc;++i) {
inlist.push_back(string(argv[i])+"/");
inlist.push_back(string(argv[i]));
if (inlist.back().rfind(".db")!=
inlist.back().length()-3) inlist.back()+="/";
else {
inlist.back()=inlist.back().substr(0,inlist.back().length()-3);
My_In_File::OpenDB(inlist.back()+"/");
}
}
MakeDir(output);
vector<string> filelist;
if (inlist[0][inlist[0].length()-1]!='/') {
sqlite3 *db=NULL;
int res=sqlite3_open((inlist[0]+".db").c_str(),&db);
if (res!=SQLITE_OK) {
msg_Error()<<METHOD<<"(): '"<<inlist[0]<<"' returns '"
<<sqlite3_errmsg(db)<<"'."<<std::endl;
}
if (db==NULL) {
msg_Error()<<METHOD<<"(): '"<<inlist[0]<<"' not found."<<std::endl;
return 1;
}
char sql[100], *zErrMsg=0;
strcpy(sql,"select file from path");
msg_IODebugging()<<METHOD<<"(\""<<inlist[0]<<"\"): {\n";
int rc=sqlite3_exec(db,sql,AddFiles,(void*)&filelist,&zErrMsg);
if(rc!=SQLITE_OK) {
msg_IODebugging()<<METHOD<<"(): '"<<inlist[0]
<<"' returns '"<<zErrMsg<<"'."<<std::endl;
sqlite3_free(zErrMsg);
sqlite3_close(db);
return 1;
}
sqlite3_close(db);
msg_IODebugging()<<"}\n";
}
else {
string flname=output+"fl.tmp";
string tmp="ls "+filter+inlist[0]+" > "+flname;
int retls(system(tmp.c_str()));
vector<string> filelist;
std::string buf;
ifstream from(flname.c_str());
while (from) {
......@@ -105,7 +148,7 @@ int main(int argc,char **argv)
}
from.close();
int retrm(system(("rm "+flname).c_str()));
}
if (check) {
double** csmatrix=new double*[inlist.size()];
int** nummatrix=new int*[inlist.size()];
......@@ -122,7 +165,7 @@ int main(int argc,char **argv)
for (i=0;i<filelist.size();i++) {
vector<Histogram*> hvec;
for (size_t j=0;j<inlist.size();j++) {
hvec.push_back(new Histogram (inlist[j]+filelist[i],1));
hvec.push_back(new Histogram (inlist[j]+"/"+filelist[i],1));
}
for (size_t j=0;j<inlist.size();j++) {
for (size_t k=j+1;k<inlist.size();k++) {
......@@ -182,12 +225,12 @@ int main(int argc,char **argv)
int sc=0;
for (i=0;i<filelist.size();i++) {
Histogram h0(inlist[0]+filelist[i],1);
Histogram h0(inlist[0]+"/"+filelist[i],1);
if (h0.Fills()>0) {
bool valid=1;
if (mode==0) h0.Restore();
for (int j=1;j<inlist.size();j++) {
Histogram histo(inlist[j]+filelist[i],1);
Histogram histo(inlist[j]+"/"+filelist[i],1);
if (histo.Fills()>0) {
if (mode==0) histo.Restore();
switch (mode) {
......@@ -209,6 +252,8 @@ int main(int argc,char **argv)
if (mode==0) h0.Finalize();
if (mode==5) h0.Scale(factor);
if (noheader) h0.SetFills(-1);
if (filelist[i].rfind("/")!=std::string::npos)
ATOOLS::MakeDir(output+filelist[i].substr(0,filelist[i].rfind("/")));
h0.Output(output+filelist[i]);
sc++;
}
......@@ -216,6 +261,10 @@ int main(int argc,char **argv)
}
if (sc>0) cout<<"successfully combined "<<inlist.size()
<<" directories containing "<<sc<<" histograms"<<endl;
for (int i(0);i<inlist.size();++i) {
if (inlist[i].rfind("/")!=inlist[i].length()-1)
My_In_File::CloseDB(inlist[i]);
}
if (rm>0) {
for (size_t i=0;i<inlist.size();++i) {
Remove(inlist[i]);
......
......@@ -4,6 +4,6 @@ bin_PROGRAMS = Combine_Analysis
Combine_Analysis_SOURCES = Combine_Analysis.C
Combine_Analysis_LDADD = @ATOOLSLIBS@ @METOOLSLIBS@ @MODELLIBS@ @PDFLIBS@ @CONDITIONAL_LHAPDFLIBS@ -lm -ldl
Combine_Analysis_LDADD = @ATOOLSLIBS@ @METOOLSLIBS@ @MODELLIBS@ @PDFLIBS@ @CONDITIONAL_LHAPDFLIBS@ @SQLITE3_LDFLAGS@ -lm -ldl
endif
......@@ -608,6 +608,7 @@ bool RootNtuple_Reader::ReadInFullEvent(Blob_List * blobs)
vars.m_type[0]=='S');
int oew(m_nlos.back()->m_n-2+(onemoreas?1:0)-vars.m_oqcd);
signalblob->SetStatus(blob_status::needs_beams);
signalblob->AddStatus(blob_status::needs_harddecays);
signalblob->SetWeight(m_weight);
// fill variation weights such that later event phases can update them
if (p_variations) {
......
......@@ -70,7 +70,23 @@ MINLO_Scale_Setter::MINLO_Scale_Setter
}
m_scale.resize(Max(m_scale.size(),m_calcs.size()));
SetCouplings();
Data_Reader read(" ",";","!","=");
Data_Reader read(" ",";","#","=");
std::vector<std::string> cores;
if (read.VectorFromFile(cores,"MINLO_ALLOW_CORE")) {
msg_Debugging()<<METHOD<<"(): Allow cores ";
Data_Reader cread(",",";","#","=");
cread.SetAddCommandLine(false);
for (size_t i(0);i<cores.size();++i) {
cread.SetString(cores[i]);
std::vector<int> flavs;
cread.VectorFromString(flavs,"");
m_cores.resize(m_cores.size()+1);
for (size_t j(0);j<flavs.size();++j)
m_cores.back().push_back(flavs[j]);
msg_Debugging()<<m_cores.back()<<" ";
}
msg_Debugging()<<"\n";
}
if (!read.ReadFromFile(m_noutmin,"MINLO_NOUT_MIN")) m_noutmin=2;
if (!read.ReadFromFile(m_cmode,"MINLO_CLUSTER_MODE")) m_cmode=1;
if (!read.ReadFromFile(m_hqmode,"MINLO_HQ_MODE")) m_hqmode=1;
......@@ -284,6 +300,31 @@ double MINLO_Scale_Setter::Calculate(const Vec4D_Vector &momenta,const size_t &m
ampl->DeleteNext();
continue;
}
if (m_cores.size()) {
bool found(true);
for (size_t i(0);i<m_cores.size();++i) {
if (ampl->Legs().size()!=m_cores[i].size()) continue;
found=false;
bool match(true);
for (size_t j(0);j<ampl->Legs().size();++j)
if (!m_cores[i][j].Includes(ampl->Leg(j)->Flav())) {
msg_Debugging()<<ampl->Leg(j)->Flav()<<" does not match "
<<m_cores[i][j]<<" in "<<m_cores[i]<<"\n";
match=false;
break;
}
if (match) {
found=true;
break;
}
}
if (!found) {
msg_Debugging()<<"core check failed\n";
ampl=ampl->Prev();
ampl->DeleteNext();
continue;
}
}
ampl->SetOrderQCD(ampl->OrderQCD()-1);
}
THROW(fatal_error,"Invalid amplitude");
......
......@@ -52,6 +52,8 @@ namespace PHASIC {
int m_vmode, m_cmode, m_hqmode, m_order, m_orderrs, m_bumode;
int m_usecomb, m_usepdfinfo, m_nlocpl, m_mufmode, m_murmode;
std::vector<ATOOLS::Flavour_Vector> m_cores;
PDF::Cluster_Param CoreScale(ATOOLS::Cluster_Amplitude *const ampl) const;
void KT2(const ATOOLS::Cluster_Amplitude *ampl,
......
......@@ -400,6 +400,7 @@ bool Initialization_Handler::InitializeTheFramework(int nr)
if (p_evtreader==NULL) THROW(fatal_error,"Event reader not found");
msg_Events()<<"SHERPA will read in the events."<<std::endl
<<" The full framework is not needed."<<std::endl;
InitializeTheHardDecays();
InitializeTheBeamRemnants();
InitializeTheIO();
InitializeTheReweighting();
......
......@@ -137,6 +137,7 @@ bool Sherpa::InitializeTheEventHandler()
if (mode==eventtype::EventReader) {
p_eventhandler->AddEventPhase(new EvtReadin_Phase(p_inithandler->GetEventReader(),
p_inithandler->GetVariations()));
p_eventhandler->AddEventPhase(new Hard_Decays(p_inithandler->GetHardDecayHandler()));
p_eventhandler->AddEventPhase(new Beam_Remnants(p_inithandler->GetBeamRemnantHandler()));
}
else {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment