From 98aba6158a1d9aa253f712f90a7b2ef3390c3db8 Mon Sep 17 00:00:00 2001 From: Jan Kotanski <kotanski@cern.ch> Date: Thu, 20 Oct 2011 07:18:47 +0200 Subject: [PATCH] adding DataModel to requirements (ReweightTools-00-00-04) --- .../GenAnalysisTools/ReweightTools/README | 0 .../ReweightTools/IPDFReweightTool.h | 34 ++ .../ReweightTools/PDFReweightTool.h | 171 ++++++ .../ReweightTools/cmt/requirements | 26 + .../share/PDFReweightTool_jobOptions.py | 44 ++ .../ReweightTools/src/PDFReweightTool.cxx | 520 ++++++++++++++++++ .../src/components/ReweightTools_entries.cxx | 12 + .../src/components/ReweightTools_load.cxx | 5 + 8 files changed, 812 insertions(+) create mode 100755 Generators/GenAnalysisTools/ReweightTools/README create mode 100644 Generators/GenAnalysisTools/ReweightTools/ReweightTools/IPDFReweightTool.h create mode 100644 Generators/GenAnalysisTools/ReweightTools/ReweightTools/PDFReweightTool.h create mode 100755 Generators/GenAnalysisTools/ReweightTools/cmt/requirements create mode 100644 Generators/GenAnalysisTools/ReweightTools/share/PDFReweightTool_jobOptions.py create mode 100644 Generators/GenAnalysisTools/ReweightTools/src/PDFReweightTool.cxx create mode 100755 Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_entries.cxx create mode 100755 Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_load.cxx diff --git a/Generators/GenAnalysisTools/ReweightTools/README b/Generators/GenAnalysisTools/ReweightTools/README new file mode 100755 index 00000000000..e69de29bb2d diff --git a/Generators/GenAnalysisTools/ReweightTools/ReweightTools/IPDFReweightTool.h b/Generators/GenAnalysisTools/ReweightTools/ReweightTools/IPDFReweightTool.h new file mode 100644 index 00000000000..7052df3fc86 --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/ReweightTools/IPDFReweightTool.h @@ -0,0 +1,34 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef IPDFREWEIGHTTOOL_H +#define IPDFREWEIGHTTOOL_H + +#include "GaudiKernel/IAlgTool.h" + +static const InterfaceID IID_IPDFReweightTool("IPDFReweightTool", 1, 0); + +/** +@class PDFReweightTool.h + +@brief Interface AlgTool base class for PDF reweighting tool + +@author Gia Khoriauli and Markus Cristinziani + Universitat Bonn + + May, 2008 +*/ +class IPDFReweightTool : virtual public IAlgTool { + public: + static const InterfaceID& interfaceID(); + +}; + + +inline const InterfaceID& IPDFReweightTool::interfaceID() +{ + return IID_IPDFReweightTool; +} + +#endif diff --git a/Generators/GenAnalysisTools/ReweightTools/ReweightTools/PDFReweightTool.h b/Generators/GenAnalysisTools/ReweightTools/ReweightTools/PDFReweightTool.h new file mode 100644 index 00000000000..9373081879b --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/ReweightTools/PDFReweightTool.h @@ -0,0 +1,171 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef PDFREWEIGHTTOOL_H +#define PDFREWEIGHTTOOL_H + +#include "GaudiKernel/AlgTool.h" +#include "ReweightTools/IPDFReweightTool.h" +#include "HepMC/GenEvent.h" +#include <string> +#include <vector> + +class StoreGateSvc; +class McEventCollection; + +//the lhapdf methods to initialize, evaluate and so on +extern "C" { + extern void initpdfset_(char*, int); + extern void initpdfsetbyname_(char*, int); + extern int numberpdf_(int&); + extern void initpdf_(int&); + extern void evolvepdf_(double&, double &, double*); + + extern void initpdfsetm_(int&,char*, int); + extern void initpdfsetbynamem_(int&, char*, int); + extern int numberpdfm_(int&,int&); + extern void initpdfm_(int&,int&); + extern void evolvepdfm_(int&,double&, double &, double*); +} + + + +/** +@class PDFReweightTool.h + +@brief tool for calculating MC event weights due to PDF uncertainties. + Can be used in systematics studies of different physics analysis + +@author Gia Khoriauli and Markus Cristinziani + Universitat Bonn + + May, 2008 +*/ +class PDFReweightTool : virtual public IPDFReweightTool, public AlgTool { + +public: + + PDFReweightTool(const std::string&, const std::string&, const IInterface*); + ~PDFReweightTool(); + + virtual StatusCode initialize(); + + //will call some lhapdf functions + //to initialize a PDF set and + //evaluate PDF values for the + //members of this PDF set + //calculates weights and + //stores them + virtual StatusCode execute(); + virtual StatusCode Reweight(HepMC::GenEvent*); + + //methods to pass interactied parton + //characteristics (flavors, x-momentum + //fractions and transmited momentum - Q) + //to this tool + void SetX1 (double); + void SetFl1 (int); + void SetX2 (double); + void SetFl2 (int); + void SetQ (double); + void SetPDFInputs(double, int, double, int, double); + + //to get them back (if necessary...) + double GetX1 (); + int GetFl1 (); + double GetX2 (); + int GetFl2 (); + double GetQ (); + + //scales m_x1 and m_x1 to the new beam energy + //if it is requested to reweight to the different c.m. energy + void ScaleXs (); + + //methods to retrive the calculated weights and their vector + double GetWeight (int); + const std::vector<double>& GetWeightVec (); + + //methods to be used by generetors to pass + //the values of evaluated PDFs to this tool + //and to retrieve those values by this tool + void SetEventPDF1 (double); + void SetEventPDF2 (double); + double GetEventPDF1 (); + double GetEventPDF2 (); + + void SetGenEvent (HepMC::GenEvent*); + HepMC::GenEvent* GetGenEvent (); + +private: + + HepMC::GenEvent* m_evt; + McEventCollection* m_mceventTESout; + StoreGateSvc* m_storeGate; + + + //propearties + bool m_GeneratorUse; + bool m_RewriteWeights; + std::string m_OrigPDFSetName; + std::string m_PDFSetName; + int m_PDFSetID; + std::string m_McEventCollection; + std::string m_McEventCollectionOut; + + bool m_DifferentCMEnergies; + double m_OrigBeamEnergy; + double m_NewBeamEnergy; + + + //evnet related variables, which are + //necesary to evaluate PDF weights + //the user of this tool has to + //initialize them in each event + int m_fl1; + int m_fl2; + double m_x1; + double m_x2; + double m_Q; + + //the vector will be filled whith + //calculated wheigts for the event + //the content of m_weights_vec + //should be stored in GenEvent + //weights vector + std::vector<double> m_weights_vec; + + //evaluated PDF values for interacting + //initial partons, which were used + //when generating event + //should also be provided by a Tool user + double m_used_pdf1; + double m_used_pdf2; + + double m_weightFactor; + +}; + + +inline void PDFReweightTool::SetGenEvent (HepMC::GenEvent* p_evt){ m_evt = p_evt; } +inline HepMC::GenEvent* PDFReweightTool::GetGenEvent () { return m_evt; } + +inline void PDFReweightTool::SetX1 (double p_x1) { m_x1 = p_x1; } +inline void PDFReweightTool::SetFl1 (int p_fl1) { m_fl1 = p_fl1; } +inline void PDFReweightTool::SetX2 (double p_x2) { m_x2 = p_x2; } +inline void PDFReweightTool::SetFl2 (int p_fl2) { m_fl2 = p_fl2; } +inline void PDFReweightTool::SetQ (double p_Q) { m_Q = p_Q; } +inline void PDFReweightTool::SetEventPDF1 (double p_used_pdf1){ m_used_pdf1 = p_used_pdf1; } +inline void PDFReweightTool::SetEventPDF2 (double p_used_pdf2){ m_used_pdf2 = p_used_pdf2; } + +inline double PDFReweightTool::GetX1 () { return m_x1; } +inline int PDFReweightTool::GetFl1 () { return m_fl1; } +inline double PDFReweightTool::GetX2 () { return m_x2; } +inline int PDFReweightTool::GetFl2 () { return m_fl2; } +inline double PDFReweightTool::GetQ () { return m_Q; } +inline double PDFReweightTool::GetEventPDF1 () { return m_used_pdf1; } +inline double PDFReweightTool::GetEventPDF2 () { return m_used_pdf2; } + +inline const std::vector<double>& PDFReweightTool::GetWeightVec () { return m_weights_vec; } + +#endif diff --git a/Generators/GenAnalysisTools/ReweightTools/cmt/requirements b/Generators/GenAnalysisTools/ReweightTools/cmt/requirements new file mode 100755 index 00000000000..a460f2c3156 --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/cmt/requirements @@ -0,0 +1,26 @@ +package ReweightTools +############################################################# +## TruthSelector: tool to select reconstructable GenParticles +############################################################# + +author Alan Poppleton <Alan.Poppleton@cern.ch> +author Sven Vahsen <SEVahsen@lbl.gov> + +use AtlasPolicy AtlasPolicy-* +use AtlasCLHEP AtlasCLHEP-* External +use GaudiInterface GaudiInterface-* External +use GeneratorObjects GeneratorObjects-* Generators +use AtlasHepMC AtlasHepMC-* External +use HepPDT v* LCG_Interfaces +use StoreGate StoreGate-* Control +use Lhapdf Lhapdf-* External +use TruthHelper TruthHelper-* Generators/GenAnalysisTools +use DataModel DataModel-* Control + +apply_pattern dual_use_library files='*.cxx' + +apply_pattern declare_joboptions files="*.py" + +private + + diff --git a/Generators/GenAnalysisTools/ReweightTools/share/PDFReweightTool_jobOptions.py b/Generators/GenAnalysisTools/ReweightTools/share/PDFReweightTool_jobOptions.py new file mode 100644 index 00000000000..0b24625919e --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/share/PDFReweightTool_jobOptions.py @@ -0,0 +1,44 @@ + +include.block("ReweightTools/PDFReweightTool_jobOptions.py") + +from AthenaCommon.AppMgr import ToolSvc +from ReweightTools.ReweightToolsConf import PDFReweightTool +ToolSvc += PDFReweightTool("pdfRwTool") + +##Default propearty values. No re-weighting will be performed +## +# If re-weigthing is during event generation then set True +ToolSvc.pdfRwTool.GeneratorUse = False +# If re-weighting to the different c.m. is required set True +ToolSvc.pdfRwTool.DifferentCMEnergies = False +# If re-qeighting between the different c.m. energies, then +# these two proppierties set the original beam energy used in +# events generation and the new beam energy to what the events +# are required to be re-weighted +ToolSvc.pdfRwTool.OrigBeamEnergy = 0. #5.0 #TeV +ToolSvc.pdfRwTool.NewBeamEnergy = 0. #3.5 #TeV +# Original PDF set name used in events generation +ToolSvc.pdfRwTool.OrigPDFSetName = '' #'cteq66.LHgrid' +# New PDF set name to which events are wanted to be re-weighted +# Leave it empty if events are required to be re-weighted using +# the original error PDF set +ToolSvc.pdfRwTool.PDFSetName = '' #MRST2006nnlo.LHgrid +# Set PDF set central value LHAGLUE number to be stored in an events +# weight container. The number of PDF set members will be stored as well +# If the new PDF set is not specified, then set the original PDF LHAGLUE value. +# Otherwise, set the new PDF set central value PDF LHAGLUE number +ToolSvc.pdfRwTool.PDFSetID = 0 #10550 #20550 +# Set the input McEventCollection name. Required if run on a POOL file +ToolSvc.pdfRwTool.McEventCollection = '' #'GEN_EVENT' #'GEN_AOD' +# Set the transient McEventCollection name. It will hold the re-weighted events for +# a superwising Athena algorithm. It can be stored in an output POOL file (if the +# supervisor algorithm produces it) or its content can be dumped into a ntuple +# Relevant when running on the POOL file. +ToolSvc.pdfRwTool.McEventCollectionOut = '' #'GEN_EVENT_PDFrw' #'GEN_AOD_PDFrw' +# Acts when re-weitghing goes on the existing McEventCollection. +# Set True if re-writing of the existing weights is required. +# typically, this is not needed and works properly for Herwig samples only +ToolSvc.pdfRwTool.RewriteWeights = False +ToolSvc.pdfRwTool.OutputLevel = ERROR + + diff --git a/Generators/GenAnalysisTools/ReweightTools/src/PDFReweightTool.cxx b/Generators/GenAnalysisTools/ReweightTools/src/PDFReweightTool.cxx new file mode 100644 index 00000000000..bf03c2d75f7 --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/src/PDFReweightTool.cxx @@ -0,0 +1,520 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +///////////////////////////////////////////////////////////// +// +// tool for calculating MC event weights due to PDF +// uncertainties. Can be used in systematics studies +// +// authors: Gia Khoriauli and Markus Cristinziani +// Universitat Bonn +// +// date: November 2007 +// +///////////////////////////////////////////////////////////// + +#include "GaudiKernel/Property.h" +#include "GaudiKernel/MsgStream.h" + +#include "StoreGate/StoreGateSvc.h" +#include "DataModel/DataVector.h" +#include "GeneratorObjects/McEventCollection.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/PdfInfo.h" +#include "ReweightTools/PDFReweightTool.h" + +#include <math.h> + +PDFReweightTool::PDFReweightTool(const std::string& type, + const std::string& name, + const IInterface* parent) + : AlgTool(type, name, parent), + m_evt(NULL), + m_mceventTESout(NULL), + m_storeGate(NULL), + m_fl1(0), + m_fl2(0), + m_x1(0.), + m_x2(0.), + m_Q(0.), + m_used_pdf1(0.), + m_used_pdf2(0.), + m_weightFactor(1.) +{ + declareProperty("GeneratorUse", m_GeneratorUse = false); + declareProperty("RewriteWeights", m_RewriteWeights = false); + declareProperty("OrigPDFSetName", m_OrigPDFSetName = ""); + declareProperty("PDFSetName", m_PDFSetName = ""); + declareProperty("PDFSetID", m_PDFSetID = 0); + declareProperty("McEventCollection", m_McEventCollection = ""); + declareProperty("McEventCollectionOut", m_McEventCollectionOut = ""); + declareProperty("DifferentCMEnergies", m_DifferentCMEnergies = false); + declareProperty("OrigBeamEnergy", m_OrigBeamEnergy = -999.); + declareProperty("NewBeamEnergy", m_NewBeamEnergy = -999.); +} + + +PDFReweightTool::~PDFReweightTool() { +} + + +StatusCode PDFReweightTool::initialize() { + MsgStream log( msgSvc(), name() ); + log << MSG::INFO + << "Initializing PDF Reweighting Tool " + << endreq; + + //retrieve StoreGate + StatusCode sc = service("StoreGateSvc", m_storeGate); + if (sc.isFailure()) { + log << MSG::ERROR + << "Unable to retrieve pointer to StoreGateSvc" + << endreq; + return sc; + } + + //consistency check for the tool propierties + //check #1 + if(m_GeneratorUse && m_DifferentCMEnergies) { + log << MSG::WARNING + <<"Reweighting between the different c.m. energies is requested during the event generation... " + <<"Usually, it is needed for the existing events. At this moment the run continues but it might fail. " + <<"Please, check the tool settings" + << endreq; + } + + //check #2 + if(m_DifferentCMEnergies) { + if (m_OrigBeamEnergy<=0. || m_NewBeamEnergy<=0.) { + log << MSG::ERROR + << "Reweighting between the different c.m. energies is requested. " + << "But the proper original and new beam energy values are not set. " + << "The corresponding propierty names are: 'OrigBeamEnergy' and 'NewBeamEnergy' " + << endreq; + + return StatusCode::FAILURE; + } + + //force do not store any PDF set id in this mode of reweighting + m_PDFSetID=0; + } + + //check #3 + if(m_OrigPDFSetName == "") { + log << MSG::ERROR + <<"The original (used or being used in the event generation) PDF set name has to be set. " + <<"Check the tool propierty 'OrigPDFSetName' " + << endreq; + + return StatusCode::FAILURE; + } + + //check #4 + if( (m_GeneratorUse || m_DifferentCMEnergies) && m_PDFSetName != "") { + log <<MSG::WARNING + <<"In case of reweighting with the event generation or the different c.m. energies modes " + <<"the second PDF set name has to be empty (it is not needed). Force 'PDFSetName' = ''. " + <<"If none of the above modes are assumed, then make the both propierties false: " + <<"'GeneratorUse' and 'DifferentCMEnergies' and re-run " + << endreq; + + m_PDFSetName = ""; + } + + //initialize PDF error sets. + int oldset=1; //do not change the number + char pdf_set_name[500]; + strcpy(pdf_set_name, m_OrigPDFSetName.c_str()); + initpdfsetbynamem_(oldset, pdf_set_name, strlen(pdf_set_name)); + log << MSG::INFO << " PDF set "<<m_OrigPDFSetName<<" has been initialized"<< endreq; + + if(m_PDFSetName != "") { + int newset=2; //do not change the number + char new_pdf_set_name[500]; + strcpy(new_pdf_set_name, m_PDFSetName.c_str()); + initpdfsetbynamem_(newset, new_pdf_set_name, strlen(new_pdf_set_name)); + log << MSG::INFO << " PDF set "<<m_PDFSetName<<" has been initialized"<< endreq; + } + + //the re-weighting 'flavor' is defined. anounce about it. + //scenario #1 + if(m_DifferentCMEnergies) { + log << MSG::INFO + <<"Reweighting from the original bean energy " + << m_OrigBeamEnergy + <<" TeV to the new energy " + << m_NewBeamEnergy + <<" TeV will be performed" + << endreq; + } + //scenario #2 + //this scenario has two possible options: + // 1) reweighting runs during the event generation + // 2) run on the already existing events + else if(m_OrigPDFSetName!="" && m_PDFSetName=="") { + log << MSG::INFO + <<"Event PDF weights will be calculated for this error PDF set - " + <<m_OrigPDFSetName + <<endreq; + } + //scenario #3 + else { // m_PDFSetName!="" + log << MSG::INFO + <<"Events generated with the original (central value) PDF - " + << m_OrigPDFSetName + <<" will be reweighted to the new error PDF set - " + <<m_PDFSetName + <<endreq; + } + + return StatusCode::SUCCESS; +} + + + +StatusCode PDFReweightTool::execute() { + + MsgStream log( msgSvc(), name() ); + log << MSG::DEBUG << " Execute PDF Reweighting Tool " << endreq; + + StatusCode sc; + + if(!m_GeneratorUse && m_McEventCollection=="" && m_McEventCollectionOut!="") { + log << MSG::WARNING + << "No input McEventCollection is specified. " + << "Can't create an output collection for reweighted events. Nothing will be done" + << endreq; + return StatusCode::SUCCESS; + } + + //runs on the existing McEventCollection and to store the re-weighted events + //in a new McEventCollection is requested + if(!m_GeneratorUse && m_McEventCollectionOut!="") { + m_mceventTESout = new McEventCollection(); + sc = m_storeGate->record( m_mceventTESout, m_McEventCollectionOut); + if( sc.isFailure() ) { + log << MSG::WARNING + << "New MC event container was not recorded in TDS" + << endreq; + return StatusCode::SUCCESS; + } + log <<MSG::DEBUG << "New MC Event Container Was Successfully Recorded" << endreq; + } + + + const McEventCollection* mceventTES; + + if(!m_GeneratorUse) { //runs on a POOL file (EVNT, AOD and etc.) + if(m_McEventCollection!="") {//input McEventCollection is specified + sc = m_storeGate->retrieve( mceventTES, m_McEventCollection); + if( sc.isFailure() || !mceventTES ) { + log << MSG::WARNING + << "No MC event container found in TDS" + << endreq; + return StatusCode::SUCCESS; + } + log <<MSG::DEBUG << "MC Event Container Was Successfully Retrieved" << endreq; + + //loop over the events in a container + McEventCollection::const_iterator iter = mceventTES->begin(); + McEventCollection::const_iterator iterEnd = mceventTES->end(); + + for (; iter!=iterEnd; iter++) { + HepMC::GenEvent* evt = new HepMC::GenEvent(**iter); + sc = this->Reweight(evt); + if( sc.isFailure() ) { + log << MSG::WARNING + << "Event PDF re-weighting failed" + << endreq; + return StatusCode::SUCCESS; + } + + if(m_McEventCollectionOut!="") { + m_mceventTESout->push_back(evt); + } + } + } + else { + log <<MSG::WARNING + << "Input MC Event Container is not specified. " + << "Reweighting won't be performed" + << endreq; + } + } + else { //runs during event generation + //assumes that the tool owner generator initializes its member + //GenEvent* pointer by setting it to the current event + HepMC::GenEvent* evt = this->GetGenEvent(); + if(evt) { + sc = this->Reweight(evt); + if( sc.isFailure() ) { + log << MSG::WARNING + << "Event PDF re-weighting failed" + << endreq; + return StatusCode::SUCCESS; + } + } + else { + log <<MSG::WARNING + << "Can't retrieve the current GenEvent. " + << "Reweighting can't be performed" + << endreq; + } + } + + return StatusCode::SUCCESS; +} + + + +//StatusCode PDFReweightTool::Reweight(const HepMC::GenEvent* evt_orig) { +StatusCode PDFReweightTool::Reweight(HepMC::GenEvent* evt) { + + MsgStream log( msgSvc(), name() ); + log << MSG::DEBUG << " Reweight PDF " << endreq; + + //safety check + if (!evt) { + log << MSG::WARNING + << " Invalid pointer to a GenEvent. Nothing will be done" + << endreq; + return StatusCode::SUCCESS; + } + + HepMC::PdfInfo* pdf_info = evt->pdf_info(); + + //safety check + if (pdf_info) { + log << MSG::DEBUG << " PdfInfo has been retrieved"<< endreq; + } + else { + log << MSG::WARNING + << " PdfInfo could not be retrieved. Nothing will be done" + << endreq; + return StatusCode::SUCCESS; + } + + //get the event info about x's, flavor's and an event scale + this->SetPDFInputs( + pdf_info->x1(), + pdf_info->id1(), + pdf_info->x2(), + pdf_info->id2(), + pdf_info->scalePDF() + ); + + //if the initial partons are gluons change their PDG values + //to the right indices of the corresponding PDF values in + //the PDF array filled by evolvepdf_() - method + int fl1 = 0; + int fl2 = 0; + if(m_fl1==21) fl1 = 0; else fl1 = m_fl1; + if(m_fl2==21) fl2 = 0; else fl2 = m_fl2; + + m_weightFactor = 1.; + + //Retrieve PDF evoluation values for the initial partons. + if (!m_GeneratorUse) { + m_used_pdf1 = pdf_info->pdf1(); + m_used_pdf2 = pdf_info->pdf2(); + } + else { + //store the fake values. they will be recalculated later + //using the central value PDF of the specified error PDF set + this->SetEventPDF1( -999. ); + this->SetEventPDF2( -999. ); + } + + + //do not change the numbers + int cvpdf=0; + int oldset=1; + + //pdf values are positive and non-zero. + //if the pdf_info is empty for some reason or tool is running during the events generation, + //then the pdf values either are not retrieved or they are not stored in the pdfinfo yet. calculate them here. + if(m_used_pdf1<=0. || m_used_pdf2<=0.) { + double f1[13]; + double f2[13]; + + initpdfm_(oldset,cvpdf); + evolvepdfm_(oldset,m_x1,m_Q,f1); + evolvepdfm_(oldset,m_x2,m_Q,f2); + + this->SetEventPDF1( f1[fl1+6] ); + this->SetEventPDF2( f2[fl2+6] ); + } + + + log << MSG::DEBUG <<" Interacted partons info : " + <<" Q="<<m_Q + <<" x1="<<m_x1 + <<" fl1="<<m_fl1 + <<" pdf1="<<m_used_pdf1 + <<" x2="<<m_x2 + <<" fl2="<<m_fl2 + <<" pdf2="<<m_used_pdf2 + << endreq; + + + //clear/reset the weight vector + m_weights_vec.clear(); + + //do not change the numbers + int n_pdf=0; + int newset=2; + + //retrieve number of error PDFs in the set + //scenario #3 + if(m_PDFSetName!="") + numberpdfm_(newset, n_pdf); + //scenario #2 + else if(!m_DifferentCMEnergies) + numberpdfm_(oldset, n_pdf); + + //scenario #1 + //needs x1 and x2 to be recalculated according to the new beam energy + if(m_DifferentCMEnergies) + this->ScaleXs(); + + //loop over error PDF set and calculate event weights. + //in case of the scenario #1 n_pdf=0 and therefore, + //only the central value PDF of the specified (original) + //error PDF set will be used + for(int i=0; i<=n_pdf; i++) { + + double f1[13]; + double f2[13]; + + //scenario #3 + if(m_PDFSetName!="") { + //init i-th member of the PDF set + initpdfm_(newset,i); + //evaluate PDFs + evolvepdfm_(newset,m_x1,m_Q,f1); + evolvepdfm_(newset,m_x2,m_Q,f2); + } + //scenario #1 and #2 + else { + initpdfm_(oldset,i); + evolvepdfm_(oldset,m_x1,m_Q,f1); + evolvepdfm_(oldset,m_x2,m_Q,f2); + } + + log << MSG::DEBUG <<" Interacted partons re-weighted info: " + <<" pdf id="<<i + <<" x1="<<m_x1<<" fl1="<<m_fl1<<" pdf1="<<f1[fl1+6] + <<" x2="<<m_x2<<" fl2="<<m_fl2<<" pdf2="<<f2[fl2+6] + <<endreq; + + + //if this is the first PDF from the PDF set (the central value PDF) + if (i==0) { + //if the used PDF values are not set, then use the PDFs evaluated + //at the current PDF set central value. + if(m_used_pdf1<=0.) { //assumes that m_used_pdf2 was not set too + this->SetEventPDF1( f1[fl1+6] ); + this->SetEventPDF2( f2[fl2+6] ); + } + //if they are there, then use them to calculate the 'weighting factor' for all new weights, + //which is equal to + // (cvPDF_new(p1) * cvPDF_new(p2)) / (cvPDF_original(p1) * cvPDF_original(p2)) + else { + if (fabs(m_used_pdf1)>0.) m_weightFactor *= f1[fl1+6]/m_used_pdf1; + if (fabs(m_used_pdf2)>0.) m_weightFactor *= f2[fl2+6]/m_used_pdf2; + + //and replace original PDF values + this->SetEventPDF1( f1[fl1+6] ); + this->SetEventPDF2( f2[fl2+6] ); + } + } + + + //calculate event weight at 'this' error PDF + double weight = m_weightFactor; + if (fabs(m_used_pdf1)>0.) weight *= f1[fl1+6]/m_used_pdf1; + if (fabs(m_used_pdf2)>0.) weight *= f2[fl2+6]/m_used_pdf2; + + log << MSG::DEBUG <<" weight = " + << weight + << endreq; + + //store the weight + m_weights_vec.push_back(weight); + } + + + //if the weights are being calculated during the event generation then nothing to re-write + if (m_GeneratorUse) m_RewriteWeights = false; + + //if re-write of the already existing PDF weights is required then save only the first weight + //which comes from a generator, flush other weights, which are supposed to be PDF weights + if(m_RewriteWeights) { + if(!(evt->weights().empty())) { + log << MSG::WARNING + <<" Event weights are going to be re-writed." + <<" This option works properly only for the HERWIG samples " + <<endreq; + + double genWeight = evt->weights()[0]; + evt->weights().clear(); + if(genWeight<99.) //so, it's probably the generator event weight but not a PDF LHAGLUE number + evt->weights().push_back(genWeight); + } + } + + //Save the unique LHAGLUE number of the PDF or the PDF set, which was used for the weights calculation. + //push that PDF LHAGLUE number and the number of the PDF set members into weight vector first + if(m_PDFSetID!=0) { + evt->weights().push_back((double)m_PDFSetID); + evt->weights().push_back((double)(n_pdf+1)); + } + + //now, store the weights into the GenEvent WeightContainer + for(unsigned int i=0; i<m_weights_vec.size(); i++) { + evt->weights().push_back(m_weights_vec[i]); + } + + //re-initialize a central value PDF for the next event + int cv_pdf = 0; + if(m_PDFSetName!="") + initpdfm_(newset, cv_pdf); + else + initpdfm_(oldset, cv_pdf); + + + return StatusCode::SUCCESS; +} + + +void PDFReweightTool::SetPDFInputs( double p_x1, + int p_fl1, + double p_x2, + int p_fl2, + double p_Q) { + SetX1(p_x1); + SetFl1(p_fl1); + SetX2(p_x2); + SetFl2(p_fl2); + SetQ(p_Q); + +} + + +double PDFReweightTool::GetWeight (int i) { + if(i>-1 && i<(int)m_weights_vec.size()) + return m_weights_vec[i]; + else + return -999.; +} + +void PDFReweightTool::ScaleXs() { + //the tool checks that m_NewBeamEnergy>0. + //in its initialize() - method + m_x1 *= m_OrigBeamEnergy/m_NewBeamEnergy; + m_x2 *= m_OrigBeamEnergy/m_NewBeamEnergy; +} + + diff --git a/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_entries.cxx b/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_entries.cxx new file mode 100755 index 00000000000..81f20c9f5c3 --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_entries.cxx @@ -0,0 +1,12 @@ +// $Id: ReweightTools_entries.cxx 349039 2011-03-03 12:48:56Z kotanski $ + +#include "GaudiKernel/DeclareFactoryEntries.h" +#include "ReweightTools/PDFReweightTool.h" + +DECLARE_TOOL_FACTORY ( PDFReweightTool ) + +DECLARE_FACTORY_ENTRIES( ReweightTools ) +{ + DECLARE_TOOL ( PDFReweightTool ) +} + diff --git a/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_load.cxx b/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_load.cxx new file mode 100755 index 00000000000..d50b9d1a430 --- /dev/null +++ b/Generators/GenAnalysisTools/ReweightTools/src/components/ReweightTools_load.cxx @@ -0,0 +1,5 @@ +// $Id: ReweightTools_load.cxx 348844 2011-03-02 15:21:05Z kotanski $ + +#include "GaudiKernel/LoadFactoryEntries.h" + +LOAD_FACTORY_ENTRIES( ReweightTools ) -- GitLab