diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/CMakeLists.txt b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3510a90caa617b8f7f3f7141b1f0ad20d81085dc
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/CMakeLists.txt
@@ -0,0 +1,59 @@
+################################################################################
+# Package: PileupReweighting
+################################################################################
+
+# Declare the package name:
+atlas_subdir( PileupReweighting )
+
+# Declare the package's dependencies:
+atlas_depends_on_subdirs( PUBLIC
+ Control/AthToolSupport/AsgTools
+ Event/xAOD/xAODEventInfo
+ GaudiKernel
+ PhysicsAnalysis/AnalysisCommon/PATInterfaces
+ PRIVATE
+ Control/AthenaBaseComps
+ Tools/PathResolver )
+
+# External dependencies:
+find_package( ROOT COMPONENTS MathCore RIO Cint Core Tree Hist pthread MathMore Minuit Minuit2 Matrix Physics HistPainter Rint )
+
+# tag ROOTBasicLibs was not recognized in automatic conversion in cmt2cmake
+
+# tag ROOTCintexLibs was not recognized in automatic conversion in cmt2cmake
+
+# Component(s) in the package:
+atlas_add_root_dictionary( PileupReweightingLib
+ PileupReweightingLibDictSource
+ ROOT_HEADERS PileupReweighting/TPileupReweighting.h Root/LinkDef.h
+ EXTERNAL_PACKAGES ROOT )
+
+atlas_add_library( PileupReweightingLib
+ Root/*.cxx
+ ${PileupReweightingLibDictSource}
+ PUBLIC_HEADERS PileupReweighting
+ PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES AsgTools xAODEventInfo GaudiKernel PATInterfaces
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps PathResolver )
+
+atlas_add_component( PileupReweighting
+ src/components/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODEventInfo GaudiKernel PATInterfaces AthenaBaseComps PathResolver PileupReweightingLib )
+
+atlas_add_dictionary( PileupReweightingDict
+ PileupReweighting/PileupReweightingDict.h
+ PileupReweighting/selection.xml
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODEventInfo GaudiKernel PATInterfaces AthenaBaseComps PathResolver PileupReweightingLib )
+
+atlas_add_executable( testPRW
+ src/testPRW.C
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODEventInfo GaudiKernel PATInterfaces AthenaBaseComps PathResolver PileupReweightingLib )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py )
+atlas_install_joboptions( share/*.py )
+atlas_install_data( share/*.root )
+
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/IPileupReweightingTool.h b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/IPileupReweightingTool.h
index b99122c814dcebdbe2a1181f7acd06b0a7c45fc2..151760209caf6e1346b0e559d46d6744dc3c0916 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/IPileupReweightingTool.h
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/IPileupReweightingTool.h
@@ -12,28 +12,86 @@
#include "TString.h"
#include "xAODEventInfo/EventInfo.h"
-namespace Root {
-class TPileupReweighting;
-}
+#include "PATInterfaces/ISystematicsTool.h"
+
+class TH1;
namespace CP {
+class TPileupReweighting;
- class IPileupReweightingTool : public virtual asg::IAsgTool {
+ class IPileupReweightingTool : public virtual asg::IAsgTool, virtual public CP::ISystematicsTool {
ASG_TOOL_INTERFACE( CP::IPileupReweightingTool )
public:
- /// Decorates with: PileupWeight , RandomRunNumber, RandomLumiBlockNumber
- virtual void apply ( const xAOD::EventInfo* eventInfo ) = 0;
+ /// Return combined pileup weight
+ virtual float getCombinedWeight( const xAOD::EventInfo& eventInfo ) = 0;
+
+ /// Same as above, but for a 'custom weight' variable
+ virtual float getCombinedWeight( const xAOD::EventInfo& eventInfo,Double_t x, Double_t y=0. ) = 0;
+
+ /// When using UnrepresentedDataAction=2, you may want to apply this additional weight to ensure sum of weights are preserved
+ virtual float getUnrepresentedDataWeight( const xAOD::EventInfo& eventInfo ) = 0;
+
+ /// Get the dataWeight used to 'unprescale' data collected from a given trigger combination. mu_dependency is recommended to be true
+ virtual float getDataWeight( const xAOD::EventInfo& eventInfo, const TString& trigger, bool mu_dependent=true ) = 0;
+
+ /// Get a random run number for this MC event, mu_dependency is recommended ... jetetmiss seem to like it muchly
+ virtual int getRandomRunNumber( const xAOD::EventInfo& eventInfo , bool mu_dependent=true) = 0;
+
+ /// Get the mu of a lumiblock ... needed to 'correct' the number in datasets
+ virtual float getLumiBlockMu( const xAOD::EventInfo& eventInfo ) = 0;
+ /// Get the integrated lumi of a lumiblock (in pb-1)
+ virtual double getLumiBlockIntegratedLumi( const xAOD::EventInfo& eventInfo ) = 0;
+
+ /// Decorates with:
+ /// MC: PileupWeight (CombinedWeight[*UnrepresentedDataWeight if action=2]), RandomRunNumber, RandomLumiBlockNumber, PRWHash
+ /// Data: corrected_averageInteractionsPerCrossing
+ /// mu_dependent says if the mu_dependency should be used for random run numbers or the data weights. You will get random run numbers of 0 for events with zero pileup weight
+ virtual StatusCode apply ( const xAOD::EventInfo& eventInfo, bool mu_dependent=true ) = 0;
+
+
+ /// return the prw hash used for fast updates of weights at the post-processing level ... see the share/makeWeightTree.C script for usage
+ virtual ULong64_t getPRWHash( const xAOD::EventInfo& eventInfo ) = 0;
+
+ // methods that go straight to the underlying pileup tool
+
+ /// Get a random lumiblock number for the given run number
+ virtual UInt_t GetRandomLumiBlockNumber(UInt_t runNumber) = 0;
+
+ /// possible alternative to using the EventBookkeepers info ... assuming you made your PRW Config file!!
+ virtual Double_t GetSumOfEventWeights(Int_t channel) = 0;
+ virtual Double_t GetNumberOfEvents(Int_t channel) = 0;
+
+ /// Get the integrated luminosity (in pb-1) between start and end run (inclusive)
+ virtual Double_t GetIntegratedLumi(UInt_t start, UInt_t end) = 0;
+ /// Total lumi (in pb-1) for a given trigger combination .. leave blank for the unprescaled
+ virtual Double_t GetIntegratedLumi(const TString& trigger="") = 0;
+ /** similar to above, but for only the given mcRunNumber/periodNumber */
+ virtual Double_t GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end) = 0;
+ /** return fraction of lumi assigned to periodNumber (or mcRunNumber) that is between start and end data run numbers*/
+ virtual Double_t GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end) = 0;
+ /** return fraction of lumi assigned to periodNumber (or mcRunNumber) with given mu, that is between start and end data run numbers*/
+ virtual Double_t GetIntegratedLumiFraction(Int_t periodNumber, Double_t mu, UInt_t start, UInt_t end) = 0;
+
+
+ /// use these methods when generating config files
+ virtual Int_t AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end) = 0;
+ virtual Int_t SetBinning(Int_t nbinsx, Double_t* xbins, Int_t nbinsy=0, Double_t* ybins=0) = 0;
+ virtual Int_t SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy=0, Double_t ylow=0, Double_t yup=0) = 0;
+ virtual Int_t SetBinning(TH1* hist) = 0;
+
+ /// Call this method once per event when in config file generating mode and you want standard mu reweighting
+ virtual int fill( const xAOD::EventInfo& eventInfo ) = 0;
+
+ /// Use this method if you want to do a generic reweighting instead
+ virtual int fill( const xAOD::EventInfo& eventInfo, Double_t x, Double_t y=0.) = 0;
+
- /// Get the dataWeight used to 'unprescale' data collected from a given trigger combination
- virtual double dataWeight ( const xAOD::EventInfo* eventInfo, TString trigger, bool mu_dependent=true ) = 0;
+ /// Get pointer to the underlying tool - expert use only: Will require #include "PileupReweighting/TPileupReweighting.h"
+ virtual CP::TPileupReweighting* expert() = 0;
- /// Get pointer to the underlying tool - expert use only
- virtual Root::TPileupReweighting* get() = 0;
- /// execute method - retrieves the eventInfo from the storegate and decorates it for you
- virtual StatusCode execute() = 0;
};
}
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/PileupReweightingTool.h b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/PileupReweightingTool.h
index 597d125cdec68ba4e3e87dd7850d29a71fa74881..7c54e8ebe20b05cd026d7397105054f0354d792f 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/PileupReweightingTool.h
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/PileupReweightingTool.h
@@ -14,6 +14,9 @@
#include "PileupReweighting/IPileupReweightingTool.h"
#include "PileupReweighting/TPileupReweighting.h"
+#include "AsgTools/ToolHandle.h"
+#include "PATInterfaces/IWeightTool.h"
+
namespace CP {
/// Implementation of the PileupReweighting tool
@@ -21,13 +24,36 @@ namespace CP {
/// @author Will Buttinger
///
- class PileupReweightingTool : public virtual IPileupReweightingTool,
- public asg::AsgTool, public Root::TPileupReweighting {
+ class PileupReweightingTool :
+ public asg::AsgTool, public virtual CP::TPileupReweighting,public virtual IPileupReweightingTool {
/// Create a proper constructor for Athena
ASG_TOOL_CLASS( PileupReweightingTool, CP::IPileupReweightingTool )
public:
+
+ /// Get a random lumiblock number for the given run number
+ virtual UInt_t GetRandomLumiBlockNumber(UInt_t runNumber) { return m_activeTool->GetRandomLumiBlockNumber(runNumber); }
+ /// Get the integrated luminosity (in pb-1) between start and end run (inclusive)
+ virtual Double_t GetIntegratedLumi(UInt_t start, UInt_t end) { return m_activeTool->GetIntegratedLumi(start, end); }
+ /// Total lumi (in pb-1) for a given trigger combination .. leave blank for the unprescaled
+ virtual Double_t GetIntegratedLumi(const TString& trigger="") { return m_activeTool->GetIntegratedLumi(trigger); }
+ /** similar to above, but for only the given mcRunNumber/periodNumber */
+ virtual Double_t GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end) { return m_activeTool->GetIntegratedLumi(periodNumber, start, end); }
+ /** return fraction of lumi assigned to periodNumber (or mcRunNumber) that is between start and end data run numbers*/
+ virtual Double_t GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end) { return m_activeTool->GetIntegratedLumiFraction(periodNumber, start, end); }
+ /** return fraction of lumi assigned to periodNumber (or mcRunNumber) with given mu, that is between start and end data run numbers*/
+ virtual Double_t GetIntegratedLumiFraction(Int_t periodNumber, Double_t mu, UInt_t start, UInt_t end) { return m_activeTool->GetIntegratedLumiFraction(periodNumber, mu, start, end); }
+ /// use these methods when generating config files
+ virtual Int_t AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end) { return m_activeTool->AddPeriod(periodNumber, start, end); }
+ virtual Int_t SetBinning(Int_t nbinsx, Double_t* xbins, Int_t nbinsy=0, Double_t* ybins=0) { return m_activeTool->SetBinning(nbinsx, xbins, nbinsy, ybins); }
+ virtual Int_t SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy=0, Double_t ylow=0, Double_t yup=0) { return m_activeTool->SetUniformBinning(nbinsx, xlow, xup, nbinsy, ylow, yup); }
+ virtual Int_t SetBinning(TH1* hist) { return m_activeTool->SetBinning(hist); }
+ /// helpful alternative to using the EventBookkeepers info ... assuming you made your PRW Config file!!
+ virtual Double_t GetSumOfEventWeights(Int_t channel) { return m_activeTool->GetSumOfEventWeights(channel); }
+ virtual Double_t GetNumberOfEvents(Int_t channel) { return m_activeTool->GetNumberOfEvents(channel) ; }
+
+
/// Constructor for standalone usage
PileupReweightingTool( const std::string& name );
@@ -37,43 +63,77 @@ namespace CP {
/// Finalize - can call the WriteToFile for us
virtual StatusCode finalize();
- /// Uses Input and Output properties to get and decorate an eventinfo object
- virtual StatusCode execute();
-
/// Decorates with:
- /// MC: PileupWeight , RandomRunNumber, RandomLumiBlockNumber
- /// Data: data weights from trigger expressions you've specified
- virtual void apply ( const xAOD::EventInfo* eventInfo );
-
- /// Get the dataWeight used to 'unprescale' data collected from a given trigger combination
- virtual double dataWeight ( const xAOD::EventInfo* eventInfo, TString trigger, bool mu_dependent=true );
-
-
- /// Get pointer to the underlying tool - expert use only
- virtual Root::TPileupReweighting* get() { return this; }
-
+ /// MC: PileupWeight (CombinedWeight[*UnrepresentedDataWeight if action=2]), RandomRunNumber, RandomLumiBlockNumber
+ /// Data: corrected_averageInteractionsPerCrossing
+ /// mu_dependent says if the mu_dependency should be used for random run numbers or the data weights. You will get random run numbers of 0 for events with zero pileup weight
+ virtual StatusCode apply ( const xAOD::EventInfo& eventInfo, bool mu_dependent=true );
+
+ /// Return combined pileup weight
+ virtual float getCombinedWeight( const xAOD::EventInfo& eventInfo );
+
+ /// Same as above, but for a 'custom weight' variable
+ virtual float getCombinedWeight( const xAOD::EventInfo& eventInfo,Double_t x, Double_t y=0. );
+
+ /// return the prw hash used for fast updates of weights at the post-processing level ... see the share/makeWeightTree.C script for usage
+ virtual ULong64_t getPRWHash( const xAOD::EventInfo& eventInfo );
+
+ /// Get the mu of a lumiblock ... needed to 'correct' the number in datasets
+ virtual float getLumiBlockMu( const xAOD::EventInfo& eventInfo );
+ /// Get the integrated lumi of a lumiblock (in pb-1)
+ virtual double getLumiBlockIntegratedLumi(const xAOD::EventInfo& eventInfo);
+
+ /// When using UnrepresentedDataAction=2, you may want to apply this additional weight to ensure sum of weights are preserved
+ virtual float getUnrepresentedDataWeight( const xAOD::EventInfo& eventInfo );
+
+ /// Get the dataWeight used to 'unprescale' data collected from a given trigger combination. mu_dependency is recommended to be true
+ virtual float getDataWeight( const xAOD::EventInfo& eventInfo, const TString& trigger, bool mu_dependent=true );
+
+ /// Get a random run number for this MC event, using mu-dependent randomization by default ... jetetmiss seem to like it muchly
+ virtual int getRandomRunNumber( const xAOD::EventInfo& eventInfo , bool mu_dependent=true);
+
+ /// Call this method once per event when in config file generating mode and you want standard mu reweighting
+ virtual int fill( const xAOD::EventInfo& eventInfo );
+
+ /// Use this method if you want to do a generic reweighting instead
+ virtual int fill( const xAOD::EventInfo& eventInfo, Double_t x, Double_t y=0.);
+
+ /// Get pointer to the underlying tool - expert use only. Will require #include "PileupReweighting/TPileupReweighting.h"
+ virtual CP::TPileupReweighting* expert() { return m_activeTool; }
+
+
+
+ /// The ISystematicsTool methods
+ bool isAffectedBySystematic( const CP::SystematicVariation& systematic ) const;
+ CP::SystematicSet affectingSystematics() const;
+ CP::SystematicSet recommendedSystematics() const;
+ CP::SystematicCode applySystematicVariation( const CP::SystematicSet& systConfig );
+
+#ifndef XAOD_STANDALONE
+ void updateHandler(Property& /*p*/);
+#endif
private:
+ std::string m_configStream;
bool m_inConfigMode;
- //Root::TPileupReweighting* m_tool;
+ CP::TPileupReweighting *m_upTool, *m_downTool; //systematic variation instances for the reweighting
+
+ double m_upVariation; double m_downVariation;
+ CP::SystematicVariation m_systUp, m_systDown;
+ CP::TPileupReweighting* m_activeTool; //defaults to this
std::vector<std::string> m_prwFiles;
std::vector<std::string> m_lumicalcFiles;
- std::vector<std::string> m_lumicalcTriggers;
-
- std::vector<std::string> m_dataWeightTriggers;
- bool m_useMuDependent;
std::string m_prefix;
- std::string m_inputKey, m_outputKey;
- int m_unrepdataAction;
int m_defaultChannel;
- double m_dataScaleFactor;
std::string m_usePeriodConfig;
+ std::map<int, bool> m_doneConfigs;
+ bool m_noWeightsMode;
- void checkInit(const xAOD::EventInfo* eventInfo);
+ ToolHandle<IWeightTool> m_weightTool;
// MN: this prevents ROOT dict generator from complaining about lack of ClassDef()
// Note: inheriting from TObject and not having ClassDef makes this class unsuitable for I/O
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/TPileupReweighting.h b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/TPileupReweighting.h
index a224bac9292e7ed146499471f5b5cfc2ef817510..161bb9c054228134579040052e7f949fdc315545 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/TPileupReweighting.h
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/TPileupReweighting.h
@@ -7,8 +7,6 @@
#ifndef __TPILEUPREWEIGHTING__
#define __TPILEUPREWEIGHTING__
-#define PILEUPREWEIGHTINGVERSION 3 //useful maybe for other packages to detect with. Each release I'll try to remember to increment this
-
/**
@class TPileupReweighting
@brief Tool to get the calculated MC pileup weight. Also does custom weights and other useful stuff.
@@ -46,7 +44,8 @@ class TTree;
class TFile;
-namespace Root {
+
+namespace CP {
class TPileupReweighting : public TNamed {
public:
@@ -58,51 +57,61 @@ namespace Root {
public:
- /** Initialize this class once before the event loop starts
- If distribution information is provided, it is assumed to be
- for the standard pileup reweighting */
- Int_t initialize(const TString dataRootFileName="",
- const TString dataRootHistName="",
- const TString mcRootFileName="",
- const TString mcRootHistName="", int runNumber=0, int channel=0);
- Int_t Initialize();
+ /** use a hardcoded period configuration */
+ Int_t UsePeriodConfig(const TString configName);
+ /** Add a histogram binning config. To modify the pileup histo binning, use "pileup" as name */
+ Int_t SetBinning(Int_t nbinsx, Double_t* xbins, Int_t nbinsy=0, Double_t* ybins=0);
+ Int_t SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy=0, Double_t ylow=0, Double_t yup=0);
+ Int_t SetBinning(TH1* hist);
- //-------------------------------------------------------------
- //Bonus method for retrieving the int lumi in each data period
- //the MuonEfficiencyCorrections package wants this
- //-------------------------------------------------------------
- std::vector<Double_t> getIntegratedLumiVector();
+ /** Set which channel should be used as a default when specific mc channel distributions cannot be found */
+ /** default channels can now be assigned for each mc run number .. mcRunNumber=-1 is the global default */
+ void SetDefaultChannel(Int_t channel, Int_t mcRunNumber=-1);
+ Int_t GetDefaultChannel(Int_t mcRunNumber=-1);
- /** total luminosity loaded and accepted by the tool (in inverse pb) */
- Double_t GetIntegratedLumi(TString trigger="");
+ /** total luminosity loaded and accepted by the tool (in inverse pb) */
+ Double_t GetIntegratedLumi(const TString& trigger="");
/** totalMC maps should also hold the total number of entries for each channel */
inline Double_t GetNumberOfEvents(Int_t channel) {
- std::map<Int_t, Double_t>& tMap = globalNumberOfEntries["pileup"];
- if(tMap.find(channel)==tMap.end()) {
+ Period* global = m_periods[-1];
+ if(!global) return 0;
+ if(global->numberOfEntries.find(channel)==global->numberOfEntries.end()) {
Error("GetNumberOfEvents", "Unknown channel: %d",channel);
return 0;
}
- return tMap[channel];
+ return global->numberOfEntries[channel];
}
inline Double_t GetSumOfEventWeights(Int_t channel) {
- std::map<Int_t, Double_t>& tMap = globalTotals["pileup"];
- if(tMap.find(channel)==tMap.end()) {
+ Period* global = m_periods[-1];
+ if(!global) return 0;
+ if(global->sumOfWeights.find(channel)==global->sumOfWeights.end()) {
Error("GetSumOfEventWeights", "Unknown channel: %d",channel);
return 0;
}
- return tMap[channel];
+ return global->sumOfWeights[channel];
}
+
+ /** Combine two period numbers. Histos are merged and the first number will be redirected to the second (the second is created if it doesn't exist) */
+ void RemapPeriod(Int_t periodNumber1, Int_t periodNumber2);
+
+
/** return fraction of lumi assigned to periodNumber (or mcRunNumber) that is between start and end data run numbers*/
Double_t GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end);
/** return fraction of lumi assigned to periodNumber (or mcRunNumber) with given mu, that is between start and end data run numbers*/
Double_t GetIntegratedLumiFraction(Int_t periodNumber, Double_t mu, UInt_t start, UInt_t end);
/** get the amount of integrated lumi between the two given runs */
- Double_t GetIntegratedLumi(UInt_t start, UInt_t end);
+ inline Double_t GetIntegratedLumi(UInt_t start, UInt_t end) { return GetIntegratedLumi(-1,start,end); }
/** similar to above, but for only the given mcRunNumber/periodNumber */
Double_t GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end);
+ /** get integrated lumi for specific run and lumiblock number .. comes from the 'unprescaled lumi', and is in pb*/
+ Double_t GetLumiBlockIntegratedLumi(Int_t runNumber, UInt_t lb);
+
+ /** get the lumiblock mu, useful for 'updating' the mu coming from data to account for new lumitags */
+ Float_t GetLumiBlockMu(Int_t runNumber, UInt_t lb);
+
//-----------------------------------------------------
//General Tool Configuring methods
//-----------------------------------------------------
@@ -110,80 +119,66 @@ namespace Root {
inline void DisableWarnings(Bool_t in) { m_SetWarnings = !in;}
/** Indicate if additional debugging information should be output */
inline void EnableDebugging(Bool_t in) { m_debugging = in;}
- /** Set which channel should be used as a default when specific mc channel distributions cannot be found */
- inline void SetDefaultChannel(Int_t channel) { m_defaultChannel=channel;}
- /** default channels can now be assigned for each mc run number */
- inline void SetDefaultChannelByMCRunNumber(Int_t channel, Int_t mcRunNumber) { m_defaultChannelsByMCRunNumber[mcRunNumber]=channel; }
- Int_t GetDefaultChannel(Int_t mcRunNumber=0);
+
/** Set how to handle configurations where some of your data has no corresponding mc to describe it
0=Default (Throw exception), 1=Subtract lumi from normalizations (i.e. discard data), 2=keep lumi and continue*/
inline void SetUnrepresentedDataAction(Int_t action, Double_t tolerance=0.05) {
- if(action<0 || action>2) {
+ if(action<0 || action>3) {
Fatal("SetUnrepresentedDataAction","Set action=%d - I'm afraid I can't let you do that, Dave",action);
throw std::runtime_error("Throwing 2001");
}
m_unrepresentedDataAction=action;
m_unrepDataTolerance=tolerance; //applicable for action=2 mode. Default is 5%
}
+ /** return the unrepresented data fraction in a given channel .. when using action=2, you will want to scale up all MC events by 1/(1-unrepFraction) to account for missed data */
+ Double_t GetUnrepresentedDataFraction(Int_t periodNumber,Int_t channel);
+ inline Float_t GetUnrepresentedDataWeight(Int_t periodNumber,Int_t channel) {
+ if(m_unrepresentedDataAction!=2) {
+ Warning("GetUnrepresentedDataWeight","You should not be applying this weight unless the UnrepresentedDataAction=2");
+ }
+ return 1./(1. - GetUnrepresentedDataFraction(periodNumber,channel));
+ }
+ //-----------------------------------------------------
+ //Methods to veto data in Action=1 mode
+ //-----------------------------------------------------
+ Bool_t IsUnrepresentedData(Int_t runNumber, Float_t x, Float_t y=0.);
+
/** Should the tool ignore period assignments encoded in config file */
inline void IgnoreConfigFilePeriods(Bool_t in) { m_ignoreFilePeriods=in; }
/** Assign an mc RunNumber to a data period */
- Int_t AddPeriod(Int_t runNumber, UInt_t start, UInt_t end);
- /** Get the PeriodNumber of the given runNumber,start and end. If not found, returns -1 */
- Int_t GetPeriodNumber(Int_t runNumber, UInt_t start, UInt_t end);
+ Int_t AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end);
/** Get the first period number with the data run number contained - assume all periods are disconnected for this to be useful */
Int_t GetFirstFoundPeriodNumber(UInt_t runNumber);
- /** Alternative method where you give it the mcRunNumber and the 'data' runNumber */
- Int_t GetPeriodNumber(Int_t mcRunNumber, UInt_t runNumber);
- /** use a hardcoded period configuration */
- Int_t UsePeriodConfig(const TString configName);
- /** Add a histogram binning config. To modify the pileup histo binning, use "pileup" as name */
- Int_t AddBinning(const TString weightName,Int_t nbinsx, Double_t* xbins);
- Int_t AddBinning(const TString weightName,Int_t nbinsx, Double_t xlow, Double_t xup);
- Int_t AddBinning(const TString weightName,Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy, Double_t ylow, Double_t yup);
- Int_t AddBinning(const TString weightName,TH1* hist);
- /** Combine two mc RunNumbers. Histos are merged and the second number will be redirected to the first */
- void MergeMCRunNumbers(Int_t runNumber1, Int_t runNumber2, Int_t newRunNumber=0);
+
+
/** Scale the LumiMetaData mu values by this factor */
- void SetDataScaleFactors(Float_t x,Float_t y=1., Float_t z=1.) { m_dataScaleFactorX=x;m_dataScaleFactorY=y;m_dataScaleFactorZ=z; }
- void SetMCScaleFactors(Float_t x,Float_t y=1.,Float_t z=1.) { m_mcScaleFactorX=x;m_mcScaleFactorY=y;m_mcScaleFactorZ=z;}
+ inline void SetDataScaleFactors(Float_t x,Float_t y=1.) { m_dataScaleFactorX=x;m_dataScaleFactorY=y; }
+ inline void SetMCScaleFactors(Float_t x,Float_t y=1.) { m_mcScaleFactorX=x;m_mcScaleFactorY=y;}
+
+
//-----------------------------------------------------
- //Methods to load config files or histograms
+ //Methods to load config files
//-----------------------------------------------------
Int_t AddConfigFile(const TString fileName);
Int_t AddLumiCalcFile(const TString fileName, const TString trigger="None");
Int_t AddMetaDataFile(const TString fileName,const TString channelBranchName="mc_channel_number");
- Int_t AddDistribution(TH1* hist, const TString weightName, Int_t runNumber, Int_t channelNumber);
- Int_t AddDistribution(TH1* hist, Int_t runNumber, Int_t channelNumber);
- Int_t AddDistributionFromFile(const TString fileName, const TString histName, const TString weightName, Int_t runNumber=0, Int_t channelNumber=0);
- Int_t AddDistributionFromFile(const TString fileName, const TString histName, Int_t runNumber=0, Int_t channelNumber=0);
+ /** Initialize this class once before the event loop starts
+ If distribution information is provided, it is assumed to be
+ for the standard pileup reweighting */
+ Int_t Initialize();
+
//-----------------------------------------------------
//Methods to get the various weights
//-----------------------------------------------------
- Float_t GetCombinedWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber,Float_t x,Float_t y=0., Float_t z=0.);
- Float_t GetPeriodWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber);
- Float_t GetPrimaryWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber,Float_t x);
- Float_t GetSecondaryWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber,Float_t x, Float_t y);
-
-
- //inline methods for above methods with the default "pileup" weight
- inline Float_t GetCombinedWeight(Int_t periodNumber, Int_t channelNumber,Float_t x, Float_t y=0., Float_t z=0.) {
- return GetCombinedWeight("pileup",periodNumber,channelNumber,x,y,z);
- }
- inline Float_t GetPeriodWeight(Int_t periodNumber, Int_t channelNumber) {
- return GetPeriodWeight("pileup",periodNumber,channelNumber);
- }
- inline Float_t GetPrimaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x) {
- return GetPrimaryWeight("pileup",periodNumber,channelNumber,x);
- }
- inline Float_t GetSecondaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x, Float_t y) {
- return GetSecondaryWeight("pileup",periodNumber,channelNumber,x,y);
- }
+ Float_t GetCombinedWeight(Int_t periodNumber, Int_t channelNumber,Float_t x,Float_t y=0.);
+ Float_t GetPeriodWeight(Int_t periodNumber, Int_t channelNumber);
+ Float_t GetPrimaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x);
+ Float_t GetSecondaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x,Float_t y);
//-----------------------------------------------------
//Methods to work with the metadata
@@ -207,17 +202,16 @@ namespace Root {
//-----------------------------------------------------
//Methods to generate config files
//-----------------------------------------------------
- Int_t Fill(const TString weightName,Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y=0., Float_t z=0.);
- Int_t Fill(Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y=0., Float_t z=0.);
+ Int_t Fill(Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y=0.);
Int_t WriteToFile(const TString filename=""); //if no name given, will use tool name
Int_t WriteToFile(TFile* outFile);
- //-----------------------------------------------------
- //Methods to veto data in Action=1 mode
- //-----------------------------------------------------
- Bool_t IsUnrepresentedData(const TString weightName, Int_t runNumber, Float_t x, Float_t y=0., Float_t z=0.);
- Bool_t IsUnrepresentedData(Int_t runNumber, Float_t x, Float_t y=0., Float_t z=0.);
+ //
+ //Methods used to do the prwTree - for fast recalculation of prw weights when processing ntuples
+ //
+ ULong64_t GetPRWHash(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y=0.);
+ Bool_t MakeWeightTree(TString channelNumbers, TString outFile, TString hashBranch="PRWHash", TString weightBranch="PileupWeight");
//-----------------------------------------------------
//RandomDataPeriod functionality stuff
@@ -226,9 +220,8 @@ namespace Root {
void SetRandomSeed(int seed) {m_random3->SetSeed(seed);}
int GetRandomSeed() {return m_random3->GetSeed();}
/** Gets a random data run number according to the integrated lumi distribution associated to this mcRunNumber */
+ /** allows use of custom reweighting config files */
UInt_t GetRandomRunNumber(Int_t mcRunNumber);
- /** Similar to above method, but allows use of custom reweighting config files .. using weightName='pileup' has same behaviour as above */
- UInt_t GetRandomRunNumber(Int_t mcRunNumber, TString weightName);
/** Get random run number according to integrated lumi distribution for the given mu value (uses the binning of the mu histogram) */
UInt_t GetRandomRunNumber(Int_t mcRunNumber, Double_t mu);
/** Get random period number from the sub-periods assigned to this run number */
@@ -241,75 +234,56 @@ namespace Root {
//Methods for PROOF cluster merging of generated histos
//-----------------------------------------------------
Int_t Merge(TCollection *coll);
- std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > >& GetInputHistograms() { return m_inputHistograms;}
+ //*std::map<Int_t,std::map<Int_t, TH1*> > & GetInputHistograms() { return m_inputHistograms;}
//-----------------------------------------------------
//Methods to inspect the input and weighting histograms
//-----------------------------------------------------
- TH1* GetInputHistogram(const TString weightName, Int_t channelNumber, Int_t runNumber) {
- if(m_inputHistograms.find(weightName)==m_inputHistograms.end()||
- m_inputHistograms[weightName].find(channelNumber)==m_inputHistograms[weightName].end()||
- m_inputHistograms[weightName][channelNumber].find(runNumber)==m_inputHistograms[weightName][channelNumber].end()) {
- return 0;
- }
- return m_inputHistograms[weightName][channelNumber][runNumber];
- }
-
- TH1D* GetPrimaryDistribution(const TString weightName, Int_t channelNumber, Int_t periodNumber) {
- if(primaryDistributions.find(weightName)==primaryDistributions.end()||
- primaryDistributions[weightName].find(channelNumber)==primaryDistributions[weightName].end()||
- primaryDistributions[weightName][channelNumber].find(periodNumber)==primaryDistributions[weightName][channelNumber].end()) {
+ TH1* GetInputHistogram(Int_t channelNumber, Int_t periodNumber) {
+ if(m_periods.find(periodNumber)==m_periods.end()||
+ m_periods[periodNumber]->inputHists.find(channelNumber)==m_periods[periodNumber]->inputHists.end()) {
return 0;
}
- return primaryDistributions[weightName][channelNumber][periodNumber];
+ return m_periods[periodNumber]->inputHists[channelNumber];
}
- TH2D* GetSecondaryDistribution(const TString weightName, Int_t channelNumber, Int_t periodNumber) {
- if(secondaryDistributions.find(weightName)==secondaryDistributions.end()||
- secondaryDistributions[weightName].find(channelNumber)==secondaryDistributions[weightName].end()||
- secondaryDistributions[weightName][channelNumber].find(periodNumber)==secondaryDistributions[weightName][channelNumber].end()) {
+ TH1D* GetPrimaryDistribution(Int_t channelNumber, Int_t periodNumber) {
+ if(m_periods.find(periodNumber)==m_periods.end()||
+ m_periods[periodNumber]->primaryHists.find(channelNumber)==m_periods[periodNumber]->primaryHists.end()) {
return 0;
}
- return secondaryDistributions[weightName][channelNumber][periodNumber];
+ return m_periods[periodNumber]->primaryHists[channelNumber];
}
TH1D* GetPrimaryTriggerDistribution(const TString trigger, Int_t periodNumber) {
- if(m_triggerPrimaryDistributions.find(trigger)==m_triggerPrimaryDistributions.end()||
- m_triggerPrimaryDistributions[trigger].find(periodNumber)==m_triggerPrimaryDistributions[trigger].end()) {
+ if(m_periods.find(periodNumber)==m_periods.end()||
+ m_periods[periodNumber]->triggerHists.find(trigger)==m_periods[periodNumber]->triggerHists.end()) {
return 0;
}
- return m_triggerPrimaryDistributions[trigger][periodNumber];
+ return m_periods[periodNumber]->triggerHists[trigger];
}
/** Method for weighting data to account for prescales and mu bias. Use by giving the tool multiple lumicalc files, one for each trigger */
Double_t GetDataWeight(Int_t runNumber, TString trigger, Double_t x);
- Double_t GetDataWeight(Int_t runNumber, TString trigger); //version without mu dependence
+ Double_t GetDataWeight(Int_t runNumber, TString trigger);//version without mu dependence
// other methods
Bool_t IsInitialized() { return m_isInitialized; }
- //-----------------------------------------------------
- //Depricated methods
- //-----------------------------------------------------
- /** Add a data distribution histogram or ttree (do this after mc) */
- void AddDataDistribution(const TString& dataRootFileName,const TString& dataRootHistName, int runNumber=0);
- /** Add a mc distribution histogram or ttree */
- void AddMCDistribution(const TString& mcRootFileName,const TString& mcRootHistName, int runNumber=0, int channel=0);
- /** Add a custom data distribution histogram or ttree (do this after mc) */
- void AddCustomDataDistribution(const TString& customRootFileName,const TString& customRootHistName, TString customName="", int runNumber=0);
- /** Add a custom mc distribution histogram or ttree */
- void AddCustomMCDistribution(const TString& customRootFileName,const TString& customRootHistName, TString customName="", int runNumber=0, int channel=0);
- Float_t getPileupWeight( float mu, int runNumber=0, int channel=0 );
-
- private:
-
- TH1* CloneEmptyHistogram(const TString weightName, Int_t runNumber, Int_t channelNumber);
+ Int_t AddDistribution(TH1* hist, Int_t runNumber, Int_t channelNumber);
+
+ protected:
+ Int_t GetNearestGoodBin(Int_t thisMCRunNumber, Int_t bin);
+ Int_t IsBadBin(Int_t thisMCRunNumber, Int_t bin);
+
+
+ TH1* CloneEmptyHistogram(Int_t runNumber, Int_t channelNumber);
/** Normalize histograms */
void normalizeHistogram(TH1* histo);
void AddDistributionTree(TTree *tree, TFile *file);
- Int_t FactorizeDistribution(TH1* hist, const TString weightName, Int_t channelNumber, Int_t periodNumber,bool includeInMCRun,bool includeInGlobal);
+ //*Int_t FactorizeDistribution(TH1* hist, const TString weightName, Int_t channelNumber, Int_t periodNumber,bool includeInMCRun,bool includeInGlobal);
void CalculatePrescaledLuminosityHistograms(const TString trigger);
@@ -317,84 +291,34 @@ namespace Root {
Bool_t m_SetWarnings;
Bool_t m_debugging;
Bool_t m_countingMode;
- Int_t m_defaultChannel;
Int_t m_unrepresentedDataAction;
- Bool_t m_isInitialized;
- Bool_t m_lumiVectorIsLoaded;
- Float_t m_dataScaleFactorX;Float_t m_dataScaleFactorY;Float_t m_dataScaleFactorZ;
- Float_t m_mcScaleFactorX;Float_t m_mcScaleFactorY;Float_t m_mcScaleFactorZ;
+ Bool_t m_isInitialized;Bool_t m_lumiVectorIsLoaded;
+ Float_t m_dataScaleFactorX;Float_t m_dataScaleFactorY;
+ Float_t m_mcScaleFactorX;Float_t m_mcScaleFactorY;
Int_t m_nextPeriodNumber;
Bool_t m_ignoreFilePeriods;
TTree *m_metadatatree;
Double_t m_unrepDataTolerance;
Bool_t m_doGlobalDataWeight; //used in GetDataWeight to flag mu-independent version
+ Int_t m_lumicalcRunNumberOffset; //used for 'faking' a lumicalc file for run2
- //-----------------------------------------------------
- //Bonus method members
- //-----------------------------------------------------
- std::vector<Double_t> m_integratedLumiVector;
+ /** map storing the lumicalc file locations - used when building DataPileupWeights */
+ std::map<TString,TString> m_lumicalcFiles;
//-----------------------------------------------------
//Shared private data members
//-----------------------------------------------------
- /** Map from mc RunNumber to (periodNumber,[start,end]) period definitions. One mc runNumber can have many disconnected periods */
- std::map<Int_t,std::map<Int_t, std::pair<UInt_t,UInt_t> > > m_periods;
- /** Quick map back from periodNumber to mcRunNumber */
- std::map<Int_t,Int_t> m_periodToMCRun;
- /** Map of empty histograms to use for custom weights. Also holds the standard pileup histogram in "pileup" */
- std::map<TString, TH1*> m_emptyHistograms;
- /** Remappings of mc run numbers. Used through the MergeMCRunNumbers method */
- std::map<Int_t,Int_t> m_mcRemappings;
- std::map<Int_t,std::map<Int_t,Int_t> > m_mcMerges;
- //-----------------------------------------------------
- //The standard pileup private members
- //-----------------------------------------------------
- /** [weightName,channelNumber] -> total (N,L,D) */
- std::map<TString, std::map<Int_t, Double_t> > globalTotals;
- /** [weightName,channelNumber,periodNumber] -> total (N_A,L_A,D_A) */
- std::map<TString, std::map<Int_t, std::map<Int_t, Double_t> > > periodTotals;
- /** [weightName,channelNumber,periodNumber] -> 1D distribution */
- std::map<TString, std::map<Int_t, std::map<Int_t, TH1D*> > > primaryDistributions;
- /** [weightName,channelNumber,periodNumber] -> 2D distribution */
- std::map<TString, std::map<Int_t, std::map<Int_t, TH2D*> > > secondaryDistributions;
-
- /** [weightName,channelNumber] -> NumberOfEntries from histograms */
- std::map<TString, std::map<Int_t, Double_t> > globalNumberOfEntries;
-
- //data is stored in channel=-1
- std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > > m_inputHistograms;
-
- //[trigger,runNumber] -> lumicalc derived histograms used to weight data to account for prescaling, based on mu
- std::map<TString, std::map<Int_t, TH1*> > m_triggerInputHistograms;
- //[trigger,periodNumber] -> assembled into the individual periods
- std::map<TString, std::map<Int_t, TH1D*> > m_triggerPrimaryDistributions;
- /** channel metadata map */
- std::map<TString, std::map<Int_t, Double_t> > m_metadata;
- //-----------------------------------------------------
- //RandomDataPeriod functionality stuff
- //numbers generated seperately for each mc period
- //-----------------------------------------------------
- TRandom3 *m_random3;
- /** Integrated lumi in each run, split in to the different mc periods: [weightName,mcRunNumber,runNumber]->lumi */
- std::map<TString, std::map<Int_t, std::map<UInt_t, Double_t> > > dataPeriodRunTotals;
- /** Integrated lumi in each run for each bin of the different mc periods: [weightName, binNumber, mcRunNumber, runNumber]->lumi*/
-#ifndef __MAKECINT__ //this map is too deep for rootcint to handle
- std::map<TString, std::map<UInt_t, std::map<Int_t, std::map<UInt_t, Double_t> > > > dataPeriodBinnedRunTotals;
-#endif
- /** [weightName,datarunnum,binnum] -> badbin flag */
- std::map<TString, std::map<Int_t, std::map<Int_t, Bool_t> > > m_badbins;
+ /** the empty histogram used for this weight... effectively holds the configuration of the binning */
+ TH1* m_emptyHistogram;
- /** mapping for default channels by mc run number - introduced for mc12a/b combined functionality */
- std::map<Int_t,Int_t> m_defaultChannelsByMCRunNumber;
- /** lumi in each lumiblock in each run loaded. used for GetRandomLumiBlockNumber */
- std::map<UInt_t, std::map<UInt_t, Double_t> > m_lumiByRunAndLbn;
- /** map storing the lumicalc file locations - used when building DataPileupWeights */
- std::map<TString,TString> m_lumicalcFiles;
+ /** channel metadata map */
+ std::map<TString, std::map<Int_t, Double_t> > m_metadata;
+
struct CompositeTrigger {
int op;
@@ -422,7 +346,62 @@ namespace Root {
CompositeTrigger* makeTrigger(TString s);
- ClassDef(TPileupReweighting,1)
+public:
+ inline void PrintPeriods() { for(auto p : m_periods) {std::cout << p.first << " -> "; p.second->print("");} }
+ struct Period {
+ Period(Int_t _id, UInt_t _start, UInt_t _end, Int_t _defaultChannel) : id(_id),start(_start),end(_end),defaultChannel(_defaultChannel) {};
+ Int_t id;
+ UInt_t start;
+ UInt_t end;
+ Int_t defaultChannel;
+ std::map<Int_t, Int_t> inputBinRedirect;
+ std::map<Int_t, Double_t> unrepData; //indexed by channel number
+ std::vector<Period*> subPeriods;
+ std::vector<UInt_t> runNumbers; //populated with runNumbers that actually had some data available
+ std::map<Int_t, TH1*> inputHists;
+ std::map<Int_t, Double_t> sumOfWeights;
+ std::map<Int_t, Int_t> numberOfEntries;
+ std::map<Int_t, TH1D*> primaryHists; //normalized histograms, indexed by channelNumber. -1 holds the data
+ std::map<Int_t, TH2D*> secondaryHists; //semi-normalized histograms, only used in 2D reweighting
+ std::map<TString, TH1D*> triggerHists; //unnormalized ... i.e. integral should be equal to the lumi! ... should really only be filled in the channel=-1 case (i.e. data)
+ bool contains(unsigned int runNumber) {
+ if(runNumber >= start && runNumber <= end) return true;
+ for(auto p : subPeriods) if(p->contains(runNumber)) return true;
+ return false;
+ };
+ void SetDefaultChannel(Int_t channel) {
+ defaultChannel=channel;
+ for(auto p : subPeriods) p->SetDefaultChannel(channel);
+ };
+ void print(const char* prefix) {
+ std::cout << prefix << id << "[" << start << "," << end << "] : ";
+ for(auto hist : inputHists) std::cout << hist.first << " , ";
+ std::cout << std::endl;
+ for(auto p : subPeriods) p->print(TString::Format(" %s",prefix).Data()); };
+ };
+ struct Run {
+ std::map<TString,TH1*> inputHists; //key is the 'trigger'
+ std::map<Int_t,Double_t> badBins;
+ Double_t lumi; //total data in run
+ std::map<UInt_t, std::pair<Double_t,Double_t> > lumiByLbn; //key=lbn, value = <lumi,mu>
+ TH1D* muDist; //mu distribution for this run
+ };
+protected:
+ std::map<Int_t, Period*> m_periods; //periods mapped by id. -1 = the "global" period (0->9999999). Uses a pointer so can easily implement remap as two entries pointing at same period
+ std::map<UInt_t, Run> m_runs; //runs mapped by runNumber
+
+ std::map<Int_t, Double_t> unrepDataByChannel; // -1 holds the total unrepData!
+
+ std::string m_prwFilesPathPrefix;
+
+ //-----------------------------------------------------
+ //RandomDataPeriod functionality stuff
+ //numbers generated seperately for each mc period
+ //-----------------------------------------------------
+ TRandom3 *m_random3;
+
+public:
+ ClassDef(TPileupReweighting,0)
}; // End: class definition
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/selection.xml b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/selection.xml
index 5e88bb132270854112fd876882173d61a925bca1..824ad0976961128cf856251006633c4ffd056535 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/selection.xml
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/PileupReweighting/selection.xml
@@ -1,4 +1,5 @@
<lcgdict>
+ <namespace name="CP"/>
<class name="CP::PileupReweightingTool" />
<class name="CP::IPileupReweightingTool" />
</lcgdict>
\ No newline at end of file
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/LinkDef.h b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/LinkDef.h
index 1b26fea9a28e053b8a04711384525341cc46059d..765abae5de05df8103fa05d4ba738ef53fcf05ba 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/LinkDef.h
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/LinkDef.h
@@ -2,11 +2,15 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
+#ifndef __PRWLINKDEF__
+#define __PRWLINKDEF__
+
+
#include <map>
#include <string>
#include "TString.h"
-#include <PileupReweighting/TPileupReweighting.h>
-#include "PathResolver/PathResolver.h"
+#include "PileupReweighting/TPileupReweighting.h"
+
#ifdef __CINT__
@@ -14,62 +18,13 @@
#pragma link off all classes;
#pragma link off all functions;
-#pragma link C++ namespace Root ;
-#pragma link C++ class Root::TPileupReweighting+ ;
-#pragma link C++ class map<int,map<int,pair<unsigned int,unsigned int > > >+ ;
-#pragma link C++ class map<int,pair<unsigned int,unsigned int > >+ ;
-#pragma link C++ class pair<unsigned int,unsigned int>+ ;
-#pragma link C++ class map<TString,TH1*>+ ;
-#pragma link C++ class map<int,map<int,int> >+ ;
-#pragma link C++ class map<TString,map<int, double> >+ ;
-#pragma link C++ class map<int,double>+ ;
-#pragma link C++ class map<TString,map<int,map<int,double> > >+;
-#pragma link C++ class map<int,map<int,double> >+;
-#pragma link C++ class map<TString,map<int, map<int,TH1D*> > >+;
-#pragma link C++ class map<int,map<int,TH1D*> >+;
-#pragma link C++ class map<int,TH1D*>+;
-#pragma link C++ class map<TString,map<int, map<int,TH2D*> > >+;
-#pragma link C++ class map<int,map<int,TH2D*> >+;
-#pragma link C++ class map<int,TH2D*>+;
-#pragma link C++ class map<TString,map<int, map<int,TH1*> > >+;
-#pragma link C++ class map<int,map<int,TH1*> >+;
-#pragma link C++ class map<int,TH1*>+;
-//#pragma link C++ class map<TString,map<unsigned int,map<int,map<unsigned int,double> > > >+;
-//#pragma link C++ class map<unsigned int,map<int,map<unsigned int,double> > >+;
-#pragma link C++ class map<TString,map<int,map<unsigned int,double> > >+;
-#pragma link C++ class map<int,map<unsigned int,double> >+;
-#pragma link C++ class map<unsigned int,double>+;
-#pragma link C++ class map<TString,map<int,map<int,bool> > >+;
-#pragma link C++ class map<int,map<int,bool> >+;
-#pragma link C++ class map<int,bool>+;
-#pragma link C++ class map<TString,TString>+;
+#pragma link C++ namespace CP;
-#pragma link C++ class pair<int,map<int,pair<unsigned int,unsigned int > > >+ ;
-#pragma link C++ class pair<int,pair<unsigned int,unsigned int > >+ ;
-#pragma link C++ class pair<TString,TH1*>+ ;
-#pragma link C++ class pair<int,map<int,int> >+ ;
-#pragma link C++ class pair<TString,map<int, double> >+ ;
-#pragma link C++ class pair<TString,map<int,map<int,double> > >+;
-#pragma link C++ class pair<int,map<int,double> >+;
-#pragma link C++ class pair<TString,map<int, map<int,TH1D*> > >+;
-#pragma link C++ class pair<int,map<int,TH1D*> >+;
-#pragma link C++ class pair<int,TH1D*>+;
-#pragma link C++ class pair<TString,map<int, map<int,TH2D*> > >+;
-#pragma link C++ class pair<int,map<int,TH2D*> >+;
-#pragma link C++ class pair<int,TH2D*>+;
-#pragma link C++ class pair<TString,map<int, map<int,TH1*> > >+;
-#pragma link C++ class pair<int,map<int,TH1*> >+;
-#pragma link C++ class pair<int,TH1*>+;
-//#pragma link C++ class pair<TString,map<unsigned int,map<int,map<unsigned int,double> > > >+;
-//#pragma link C++ class pair<unsigned int,map<int,map<unsigned int,double> > >+;
-#pragma link C++ class pair<TString,map<int,map<unsigned int,double> > >+;
-#pragma link C++ class pair<int,map<unsigned int,double> >+;
-#pragma link C++ class pair<unsigned int,double>+;
-#pragma link C++ class pair<TString,map<int,map<int,bool> > >+;
-#pragma link C++ class pair<int,map<int,bool> >+;
-#pragma link C++ class pair<int,bool>+;
+#pragma link C++ class CP::TPileupReweighting+ ;
+#endif
+
#endif
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/PileupReweightingTool.cxx b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/PileupReweightingTool.cxx
index dae691f71a9660904e9060391880890f22475b04..0accfd4d1f41e9860b5ea0135d2c78ede3133e7a 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/PileupReweightingTool.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/PileupReweightingTool.cxx
@@ -9,35 +9,144 @@
#include "PathResolver/PathResolver.h"
-#include "xAODEventInfo/EventInfoAuxContainer.h"
+#include "PATInterfaces/SystematicRegistry.h"
+
+#ifndef XAOD_STANDALONE
+#include "GaudiKernel/ITHistSvc.h"
+#endif
+
+#ifdef XAOD_STANDALONE
+#include "ReweightUtils/McEventWeight.h"
+#endif
+
+#include "TH1.h"
+#include "TTree.h"
+
namespace CP {
-PileupReweightingTool::PileupReweightingTool( const std::string& name ) : asg::AsgTool( name ), Root::TPileupReweighting(name.c_str()), m_inConfigMode(false) {
+PileupReweightingTool::PileupReweightingTool( const std::string& name ) :CP::TPileupReweighting(name.c_str()), asg::AsgTool( name ),
+ m_inConfigMode(false),
+ m_upTool(0), m_downTool(0), m_systUp("PRW_DATASF", 1 ), m_systDown("PRW_DATASF", -1),
+ m_activeTool(this),
+ m_noWeightsMode(false),
+ m_weightTool("McEventWeight/myWeightTool") {
+
+#ifndef XAOD_STANDALONE
+ declareProperty("ConfigOutputStream", m_configStream="", "When creating PRW config files, this is the THistSvc stream it goes into. If blank, it wont write this way");
+#endif
declareProperty("ConfigFiles", m_prwFiles, "List of prw config files"); //array of files
- declareProperty("LumiCalcFiles", m_lumicalcFiles, "List of lumicalc files"); //array of files
- declareProperty("LumiCalcTriggers",m_lumicalcTriggers, "List of triggers associated to each lumicalc file. 'None' is a keyword. If none specified, will assume all are 'None'"); //alias names - if not same length of LumiCalcFiles, will just use "None" (the default)
- declareProperty("UnrepresentedDataAction",m_unrepdataAction=2);
- declareProperty("DataScaleFactor",m_dataScaleFactor=1.);
+ declareProperty("ConfigFilesPathPrefix", m_prwFilesPathPrefix="", "Path of additional folder structure in prw config files"); //string prefix
+ declareProperty("LumiCalcFiles", m_lumicalcFiles, "List of lumicalc files, in the format '<filename>:<trigger>' .. if no trigger given, 'None' is assumed"); //array of files
+ declareProperty("Prefix",m_prefix="","Prefix to attach to all decorations ... only used in the 'apply' method");
+ declareProperty("UnrepresentedDataAction",m_unrepresentedDataAction=3,"1 = remove unrepresented data, 2 = leave it there, 3 = reassign it to nearest represented bin");
+ declareProperty("DataScaleFactor",m_dataScaleFactorX=1.);
declareProperty("DefaultChannel",m_defaultChannel=0); //when channel info not present in config file, use this channel instead
- declareProperty("UsePeriodConfig",m_usePeriodConfig="","Use this period configuration when in config generating mode. Leave blank to try auto-detect");
- declareProperty("DataWeightTriggers",m_dataWeightTriggers,"List of triggers used to calculate dataweights and attach to eventinfo during apply method, if running on data");
- declareProperty("MuDependentDataWeights",m_useMuDependent=true,"Mu dependent data weights means calculating dataweights from luminosity in the bins of mu - more accurate than doing a global weight with all the luminosity");
- declareProperty("Prefix",m_prefix="","Prefix to attach to all decorations");
- declareProperty("Input",m_inputKey="","Specify a specific EventInfo object");
- declareProperty("Output",m_outputKey="","Specify an output EventInfo object. If differs from input, will create a clone of EventInfo and decorate that");
+ declareProperty("UsePeriodConfig",m_usePeriodConfig="auto","Use this period configuration when in config generating mode. Set to 'auto' to auto-detect");
+
+ declareProperty("DataScaleFactorUP",m_upVariation=0.,"Set to a value representing the 'up' fluctuation - will report a PRW_DATASF uncertainty to Systematic Registry");
+ declareProperty("DataScaleFactorDOWN",m_downVariation=0.,"Set to a value representing the 'down' fluctuation - will report a PRW_DATASF uncertainty to Systematic Registry");
+
+ declareProperty("LumiCalcRunNumberOffset",m_lumicalcRunNumberOffset=0,"Use to 'fake' a Run2 lumicalc file. Suggest using a value of 22000 to do this from an 8TeV lumicalc file");
+
+#ifdef XAOD_STANDALONE
+ declareProperty( "WeightTool", m_weightTool = new McEventWeight("myWeightTool"),"The tool to compute the weight in the sumOfWeights");
+#else
+ declareProperty( "WeightTool", m_weightTool,"The tool to compute the weight in the sumOfWeights");
+#endif
+
+#ifndef XAOD_STANDALONE
+ //attached update handler to the outputlevel property, so we can pass changes on to the underlying tool
+ auto props = getProperties();
+ for( Property* prop : props ) {
+ if( prop->name() != "OutputLevel" ) {
+ continue;
+ }
+ prop->declareUpdateHandler( &PileupReweightingTool::updateHandler, this );
+ break;
+ }
+#endif
+
+}
+
+#ifndef XAOD_STANDALONE
+//rootcore can't do this yet!
+void PileupReweightingTool::updateHandler(Property& p) {
+ //call the original update handler
+ this->msg_update_handler(p);
+ EnableDebugging(this->msgLvl(MSG::DEBUG));
+}
+#endif
+
+bool PileupReweightingTool::isAffectedBySystematic( const CP::SystematicVariation& systematic ) const {
+ CP::SystematicSet sys = affectingSystematics(); return sys.find( systematic ) != sys.end();
+}
+
+CP::SystematicSet PileupReweightingTool::affectingSystematics() const {
+ CP::SystematicSet result;
+ if(m_upVariation) result.insert( m_systUp );
+ if(m_downVariation) result.insert( m_systDown );
+ return result;
+}
+
+CP::SystematicSet PileupReweightingTool::recommendedSystematics() const {
+ return affectingSystematics();
+}
+
+CP::SystematicCode PileupReweightingTool::applySystematicVariation( const CP::SystematicSet& systConfig ) {
+ if(systConfig.find( m_systUp ) != systConfig.end() && systConfig.find( m_systDown ) != systConfig.end()) {
+ ATH_MSG_ERROR("Errr... what the!? You're trying to do both PRW_DATASF directions at the same time!!!??");
+ return SystematicCode::Unsupported;
+ }
+ if(systConfig.find( m_systUp ) != systConfig.end()) {
+ if(!m_upTool) { ATH_MSG_ERROR("Requested up variation of PRW_DATASF, but not configured to do this :-("); return SystematicCode::Unsupported; }
+ m_activeTool = m_upTool;
+ }
+ else if(systConfig.find( m_systDown ) != systConfig.end() ) {
+ if(!m_downTool) { ATH_MSG_ERROR("Requested down variation of PRW_DATASF, but not configured to do this :-("); return SystematicCode::Unsupported; }
+ m_activeTool = m_downTool;
+ }
+ else m_activeTool = this;
+ return SystematicCode::Ok;
}
+
+float PileupReweightingTool::getLumiBlockMu( const xAOD::EventInfo& eventInfo ) { return CP::TPileupReweighting::GetLumiBlockMu(eventInfo.runNumber(),eventInfo.lumiBlock()); }
+
+double PileupReweightingTool::getLumiBlockIntegratedLumi(const xAOD::EventInfo& eventInfo) { return CP::TPileupReweighting::GetLumiBlockIntegratedLumi(eventInfo.runNumber(),eventInfo.lumiBlock()); }
+
StatusCode PileupReweightingTool::initialize() {
+ //set debugging if debugging is on:
+ EnableDebugging(this->msgLvl(MSG::DEBUG));
+
+ //see if we need variations
+ if(m_upVariation) {
+ m_upTool = new TPileupReweighting((name()+"_upVariation").c_str());
+ m_upTool->SetUnrepresentedDataAction(m_unrepresentedDataAction);
+ m_upTool->SetDataScaleFactors(m_upVariation);
+ }
+ if(m_downVariation) {
+ m_downTool = new TPileupReweighting((name()+"_downVariation").c_str());
+ m_downTool->SetUnrepresentedDataAction(m_unrepresentedDataAction);
+ m_downTool->SetDataScaleFactors(m_downVariation);
+ }
+
+ SetDefaultChannel(m_defaultChannel); if(m_upTool) m_upTool->SetDefaultChannel(m_defaultChannel); if(m_downTool) m_downTool->SetDefaultChannel(m_defaultChannel);
+
//should we set the period config (file maker mode)
if(m_prwFiles.size()+m_lumicalcFiles.size()==0) {
m_inConfigMode=true;
ATH_MSG_INFO("In Config file making mode.");
- if(m_usePeriodConfig!="") {
+ if(m_usePeriodConfig!="auto" && m_usePeriodConfig!="") {
/*m_tool->*/UsePeriodConfig(m_usePeriodConfig);
- }
+ }
+#ifndef XAOD_STANDALONE
+ //retrieve histsvc now, otherwise it wont get correctly configured by the time it's used in the finalize ... so it seems
+ ServiceHandle<ITHistSvc> histSvc("THistSvc",name());
+ CHECK( histSvc.retrieve() );
+#endif
} else {
//have we any prw to load
for(unsigned int j=0;j<m_prwFiles.size();j++) {
@@ -45,117 +154,243 @@ StatusCode PileupReweightingTool::initialize() {
std::string file = PathResolverFindCalibFile(m_prwFiles[j]);
if(file=="") { ATH_MSG_ERROR("Unable to find the PRW Config file: " << m_prwFiles[j]); return StatusCode::FAILURE; }
/*m_tool->*/AddConfigFile(file.c_str());
+ if(m_upTool) m_upTool->AddConfigFile(file.c_str()); if(m_downTool) m_downTool->AddConfigFile(file.c_str());
}
- //have we any lumicalc files to load
- bool useNoneLabel = (m_lumicalcFiles.size() != m_lumicalcTriggers.size());
+ //have we any lumicalc files to load? .. if we do and had no prwFiles then the user must specify the period configuration
+ if(m_lumicalcFiles.size()>0 && m_prwFiles.size()==0) {
+ if(m_usePeriodConfig!="auto" && m_usePeriodConfig!="") {
+ ATH_MSG_INFO("Setting up without a PRW config file, but with period config " << m_usePeriodConfig << ". You will only be able to use random run number and data weight functionality... no reweighting!");
+ if(UsePeriodConfig(m_usePeriodConfig)!=0) {
+ ATH_MSG_ERROR("Unrecognised PeriodConfig: " << m_usePeriodConfig); return StatusCode::FAILURE;
+ }
+ m_noWeightsMode=true; //will stop the prw weight being decorated in apply method
+ if(m_upTool) m_upTool->UsePeriodConfig(m_usePeriodConfig); if(m_downTool) m_downTool->UsePeriodConfig(m_usePeriodConfig);
+ } else {
+ ATH_MSG_INFO("No config files provided, but " << m_lumicalcFiles.size() << " lumicalc file provided. Please specify a UsePeriodConfig if you want to use the tool without a config file (e.g. do 'MC15') ");
+ return StatusCode::FAILURE;
+ }
+ }
+
for(unsigned int j=0;j<m_lumicalcFiles.size();j++) {
- ATH_MSG_INFO("Locating File: " << m_lumicalcFiles[j]);
- std::string file = PathResolverFindCalibFile(m_lumicalcFiles[j]);
- if(file=="") { ATH_MSG_ERROR("Unable to find the Lumicalc file: " << m_lumicalcFiles[j]); return StatusCode::FAILURE; }
- /*m_tool->*/SetDataScaleFactors(m_dataScaleFactor);
- /*m_tool->*/AddLumiCalcFile(file.c_str(),(useNoneLabel)? "None" : m_lumicalcTriggers[j]);
+ //see if there's a trigger at the end of the filename .. format is "file:trigger"
+ TString myFile = m_lumicalcFiles[j];
+ TString myTrigger = (myFile.Contains(':')) ? TString(myFile(myFile.Last(':')+1,myFile.Length()-myFile.Last(':'))) : TString("None");
+ myFile = (myFile.Contains(':')) ? TString(myFile(0,myFile.Last(':'))) : myFile;
+ ATH_MSG_INFO("Locating File: " << myFile);
+ std::string file = PathResolverFindCalibFile(myFile.Data());
+ if(file=="") { ATH_MSG_ERROR("Unable to find the Lumicalc file: " << myFile); return StatusCode::FAILURE; }
+ /*m_tool->*/AddLumiCalcFile(file.c_str(),myTrigger);
+ if(m_upTool) m_upTool->AddLumiCalcFile(file.c_str(),myTrigger); if(m_downTool) m_downTool->AddLumiCalcFile(file.c_str(),myTrigger);
}
- /*m_tool->*/SetUnrepresentedDataAction(m_unrepdataAction);
- /*m_tool->*/SetDefaultChannel(m_defaultChannel);
}
+ //register ourselves with the systematic registry!
+ CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance();
+ if( registry.registerSystematics( *this ) != CP::SystematicCode::Ok ) return StatusCode::FAILURE;
+
+
//delay initializing underlying tool until first usage, just in case user wants to do any advanced initialization options
+ ATH_MSG_DEBUG("Retrieving weight tool...");
+ StatusCode sc = m_weightTool.retrieve();
+ if (sc != StatusCode::SUCCESS) {
+ ATH_MSG_ERROR(" " << m_weightTool->name() << " could not be retrieved.");
+ return sc;
+ }
+ else {
+ ATH_MSG_DEBUG(" " << m_weightTool->name() << " retrieved.");
+ m_weightTool->print();
+ }
+
return StatusCode::SUCCESS;
}
StatusCode PileupReweightingTool::finalize() {
if(m_inConfigMode) {
+#ifndef XAOD_STANDALONE
+ if(m_configStream=="") {
//write the prw config files
std::string nameWithoutParent = this->name().substr(this->name().find(".")+1);
TString fileName = TString::Format("%s.prw.root", nameWithoutParent.c_str());
/*m_tool->*/WriteToFile(fileName);
- }
- return StatusCode::SUCCESS;
-}
+ } else {
+ //write to the histsvc stream instead ...
+
+ ServiceHandle<ITHistSvc> histSvc("THistSvc",name());
+ CHECK( histSvc.retrieve() );
+
+ TTree *outTreeMC=0;
+ TTree *outTreeData=0;
+ Int_t channel = 0;UInt_t runNumber = 0;
+ std::vector<UInt_t>* pStarts = 0;std::vector<UInt_t>* pEnds = 0;
+ Char_t histName[150];
+
+
+ //loop over periods ... periods only get entry in table if they have an input histogram
+ for(auto period : m_periods) {
+ if(!period.second) continue; //should never happen, but just in case!
+ if(period.first != period.second->id) continue; //skips redirects
+ runNumber = period.first;
+ pStarts = new std::vector<UInt_t>;pEnds = new std::vector<UInt_t>;
+ if(period.second->subPeriods.size()==0) {
+ pStarts->push_back(period.second->start); pEnds->push_back(period.second->end);
+ }
+ else {
+ for(auto subp : period.second->subPeriods) {
+ pStarts->push_back(subp->start); pEnds->push_back(subp->end);
+ }
+ }
+ for(auto inHist : period.second->inputHists) {
+ channel = inHist.first;
+ TH1* hist = inHist.second;
+ strcpy(histName,hist->GetName());
+ CHECK( histSvc->regHist(TString::Format("/%s/PileupReweighting/%s",m_configStream.c_str(),hist->GetName()).Data(),hist) );
+ if(!outTreeMC) {
+ outTreeMC = new TTree("MCPileupReweighting","MCPileupReweighting");
+ outTreeMC->Branch("Channel",&channel);
+ outTreeMC->Branch("RunNumber",&runNumber);
+ outTreeMC->Branch("PeriodStarts",&pStarts);
+ outTreeMC->Branch("PeriodEnds",&pEnds);
+ outTreeMC->Branch("HistName",&histName,"HistName/C");
+ CHECK( histSvc->regTree(TString::Format("/%s/PileupReweighting/%s",m_configStream.c_str(),outTreeMC->GetName()).Data(),outTreeMC) );
+ }
+ outTreeMC->Fill();
+ }
+ delete pStarts; delete pEnds;
+ }
-void PileupReweightingTool::checkInit(const xAOD::EventInfo* eventInfo) {
- if(!/*m_tool->*/IsInitialized()) {
- if(m_usePeriodConfig=="") { //try autodetect based on runnum of first event
- switch(eventInfo->runNumber()) {
- case 180164: case 183003: case 186169: case 189751: m_usePeriodConfig="MC12c";/*m_tool->*/UsePeriodConfig("MC12c");break;
- case 195847: case 195848: m_usePeriodConfig="MC12ab";/*m_tool->*/UsePeriodConfig("MC12ab");break;
- case 212272: m_usePeriodConfig="MC14";/*m_tool->*/UsePeriodConfig("MC14");break;
+ //loop over data
+ for(auto& run : m_runs) {
+ runNumber = run.first;
+ if(run.second.inputHists.find("None")==run.second.inputHists.end()) continue;
+
+ TH1* hist = run.second.inputHists["None"];
+ strcpy(histName,hist->GetName());
+ CHECK( histSvc->regHist(TString::Format("/%s/PileupReweighting/%s",m_configStream.c_str(),hist->GetName()).Data(),hist) );
+ if(!outTreeData) {
+ outTreeData = new TTree("DataPileupReweighting","DataPileupReweighting");
+ outTreeData->Branch("RunNumber",&runNumber);
+ outTreeData->Branch("HistName",&histName,"HistName/C");
+ CHECK( histSvc->regTree(TString::Format("/%s/PileupReweighting/%s",m_configStream.c_str(),outTreeData->GetName()).Data(),outTreeData) );
}
+ outTreeData->Fill();
}
- if(/*m_tool->*/Initialize()!=0) throw std::runtime_error("PileupReweighting: Failure in init");
- }
- return;
-}
-void PileupReweightingTool::apply(const xAOD::EventInfo* eventInfo) {
- checkInit(eventInfo);
+ Info("WriteToFile", "Successfully generated config file to stream: %s",m_configStream.c_str());
+ Info("WriteToFile", "Happy Reweighting :-)");
- if(!eventInfo->eventType(xAOD::EventInfo::IS_SIMULATION)) {
- //allowed to do data weights
- for(auto& trigger : m_dataWeightTriggers) {
- eventInfo->auxdecor<double>(m_prefix+trigger) = dataWeight(eventInfo, trigger, m_useMuDependent);
- }
- return;
}
- if(m_inConfigMode) {
- Fill(eventInfo->runNumber(), eventInfo->mcChannelNumber(), eventInfo->mcEventWeight(), eventInfo->averageInteractionsPerCrossing());
- return;
+#else
+ //write the prw config files
+ std::string nameWithoutParent = this->name().substr(this->name().find(".")+1);
+ TString fileName = TString::Format("%s.prw.root", nameWithoutParent.c_str());
+ /*m_tool->*/WriteToFile(fileName);
+#endif
}
- //decorate with standard PileupWeight
- eventInfo->auxdecor<double>(m_prefix+"PileupWeight") = /*m_tool->*/GetCombinedWeight( eventInfo->runNumber(), eventInfo->mcChannelNumber(), eventInfo->averageInteractionsPerCrossing() );
- //decorate with random run number etc
- //using mu-dependent randomization by default ... jetetmiss seem to like it muchly
- /*m_tool->*/SetRandomSeed(314159+eventInfo->mcChannelNumber()*2718+eventInfo->eventNumber()); //to make randomization repeatable
- eventInfo->auxdecor<unsigned int>(m_prefix+"RandomRunNumber") = /*m_tool->*/GetRandomRunNumber( eventInfo->runNumber(), eventInfo->averageInteractionsPerCrossing() );
- eventInfo->auxdecor<unsigned int>(m_prefix+"RandomLumiBlockNumber") = (eventInfo->auxdecor<unsigned int>(m_prefix+"RandomRunNumber")==0) ? 0 : /*m_tool->*/GetRandomLumiBlockNumber( eventInfo->auxdecor<unsigned int>(m_prefix+"RandomRunNumber") );
+ if(m_upTool) {delete m_upTool;m_upTool=0;} if(m_downTool) {delete m_downTool;m_downTool=0;}
- return; //FIXME return statuscode instead maybe...
+ return StatusCode::SUCCESS;
}
-double PileupReweightingTool::dataWeight(const xAOD::EventInfo* eventInfo, TString trigger, bool mu_dependent) {
- checkInit(eventInfo);
+ULong64_t PileupReweightingTool::getPRWHash( const xAOD::EventInfo& eventInfo ) {
+ return m_activeTool->GetPRWHash(eventInfo.runNumber(), eventInfo.mcChannelNumber(), eventInfo.averageInteractionsPerCrossing() );
+}
- if(eventInfo->eventType(xAOD::EventInfo::IS_SIMULATION)) {
- ATH_MSG_WARNING("Requesting Data Weight is not intended for simulated events. Returning 1.");
- return 1; //no data weights for simulated events
+/// Return combined pileup weight
+float PileupReweightingTool::getCombinedWeight( const xAOD::EventInfo& eventInfo ) {
+ return m_activeTool->GetCombinedWeight( eventInfo.runNumber(), eventInfo.mcChannelNumber(), eventInfo.averageInteractionsPerCrossing() );
+}
+
+float PileupReweightingTool::getCombinedWeight( const xAOD::EventInfo& eventInfo, Double_t x, Double_t y ) {
+ return m_activeTool->GetCombinedWeight( eventInfo.runNumber(), eventInfo.mcChannelNumber(), x, y );
+}
+
+/// When using UnrepresentedDataAction=2, you may want to apply this additional weight to ensure sum of weights are preserved
+float PileupReweightingTool::getUnrepresentedDataWeight( const xAOD::EventInfo& eventInfo ) {
+ return m_activeTool->GetUnrepresentedDataWeight( eventInfo.runNumber(), eventInfo.mcChannelNumber() );
+}
+
+
+/// Get a random run number for this MC event, mu_dependency is not recommended for this
+int PileupReweightingTool::getRandomRunNumber( const xAOD::EventInfo& eventInfo , bool mu_dependent) {
+ m_activeTool->SetRandomSeed(314159+eventInfo.mcChannelNumber()*2718+eventInfo.eventNumber()); //to make randomization repeatable
+ if(mu_dependent) return m_activeTool->GetRandomRunNumber( eventInfo.runNumber(), eventInfo.averageInteractionsPerCrossing() );
+ else return m_activeTool->GetRandomRunNumber( eventInfo.runNumber() );
+}
+
+int PileupReweightingTool::fill( const xAOD::EventInfo& eventInfo ) {
+ return fill(eventInfo,eventInfo.averageInteractionsPerCrossing());
+}
+
+int PileupReweightingTool::fill( const xAOD::EventInfo& eventInfo, Double_t x, Double_t y) {
+
+ //auto-detect the period config if necessary
+ if(m_usePeriodConfig=="auto" && !m_doneConfigs[eventInfo.runNumber()]) { //try autodetect based on runnum of first event ...
+ //if data, we only need to ensure a binning is done ... for now, we assume the MC15 binning
+ if(!eventInfo.eventType(xAOD::EventInfo::IS_SIMULATION)) {
+ if(!m_emptyHistogram) { SetUniformBinning(100,0,100.); } //use a default binning
+ } else {
+ switch(eventInfo.runNumber()) {
+ //case 180164: case 183003: case 186169: case 189751: m_usePeriodConfig="MC12c";/*m_tool->*/UsePeriodConfig("MC12c");break;
+ //case 195847: case 195848: m_usePeriodConfig="MC12ab";/*m_tool->*/UsePeriodConfig("MC12ab");break;
+ case 212272: UsePeriodConfig("MC14_8TeV");break;
+ case 222222: UsePeriodConfig("MC14_13TeV");break;
+ case 222510: case 222525: UsePeriodConfig("MC15"); break;
+ }
+ }
+ m_doneConfigs[eventInfo.runNumber()] = true;
}
- if(!mu_dependent) return /*m_tool->*/GetDataWeight(eventInfo->runNumber(), trigger);
- return /*m_tool->*/GetDataWeight( eventInfo->runNumber(), trigger, eventInfo->averageInteractionsPerCrossing() );
+ return TPileupReweighting::Fill(eventInfo.runNumber(), eventInfo.eventType(xAOD::EventInfo::IS_SIMULATION) ? eventInfo.mcChannelNumber() : -1 /*data*/, eventInfo.eventType(xAOD::EventInfo::IS_SIMULATION) ? m_weightTool->getWeight() : 1., x, y);
}
-StatusCode PileupReweightingTool::execute() {
- const xAOD::EventInfo* evtInfo =0;
+StatusCode PileupReweightingTool::apply(const xAOD::EventInfo& eventInfo, bool mu_dependent) {
- if(m_inputKey.length()>0) ATH_CHECK(evtStore()->retrieve(evtInfo,m_inputKey));
- else {
- #ifdef XAOD_STANDALONE
- ATH_CHECK( evtStore()->retrieve(evtInfo,"") ); //apparently TEvent can't do keyless retrieves!!??
- #else
- ATH_CHECK(evtStore()->retrieve(evtInfo));
- #endif
+ if(m_inConfigMode) {
+ fill( eventInfo );
+ return StatusCode::SUCCESS;
}
- //do we need to make a copy??
- if(m_inputKey!=m_outputKey && m_outputKey!="") {
- xAOD::EventInfo* evtInfoCopy = new xAOD::EventInfo( *evtInfo );
- xAOD::EventInfoAuxContainer* aux = new xAOD::EventInfoAuxContainer;
- evtInfoCopy->setStore(aux);
- ATH_CHECK( evtStore()->record( evtInfoCopy , m_outputKey ) );
- ATH_CHECK( evtStore()->record( aux , m_outputKey+"Aux." ) );
- evtInfo = evtInfoCopy;
+ if(!eventInfo.eventType(xAOD::EventInfo::IS_SIMULATION)) {
+ eventInfo.auxdecor<float>(m_prefix+"corrected_averageInteractionsPerCrossing") = getLumiBlockMu(eventInfo);
+ return StatusCode::SUCCESS;
}
- apply(evtInfo);
+ //just copy the value over for MC
+ eventInfo.auxdecor<float>(m_prefix+"corrected_averageInteractionsPerCrossing") = eventInfo.averageInteractionsPerCrossing();
+
+ //decorate with standard PileupWeight
+ if(!m_noWeightsMode) {
+ double weight = /*m_tool->*/getCombinedWeight( eventInfo );
+ if(m_unrepresentedDataAction==2) weight *= getUnrepresentedDataWeight( eventInfo );
+ eventInfo.auxdecor<double>(m_prefix+"PileupWeight") = weight;
+ }
+ //decorate with random run number etc
+ eventInfo.auxdecor<unsigned int>(m_prefix+"RandomRunNumber") = /*m_tool->*/getRandomRunNumber( eventInfo, mu_dependent );
+ eventInfo.auxdecor<unsigned int>(m_prefix+"RandomLumiBlockNumber") = (eventInfo.auxdecor<unsigned int>(m_prefix+"RandomRunNumber")==0) ? 0 : /*m_tool->*/GetRandomLumiBlockNumber( eventInfo.auxdecor<unsigned int>(m_prefix+"RandomRunNumber") );
+
+ eventInfo.auxdecor<ULong64_t>(m_prefix+"PRWHash") = getPRWHash( eventInfo );
+
+ ATH_MSG_VERBOSE("PileupWeight = " << eventInfo.auxdecor<double>(m_prefix+"PileupWeight") << " RandomRunNumber = " << eventInfo.auxdecor<unsigned int>(m_prefix+"RandomRunNumber") << " RandomLumiBlockNumber = " << eventInfo.auxdecor<unsigned int>(m_prefix+"RandomLumiBlockNumber"));
return StatusCode::SUCCESS;
+}
+
+float PileupReweightingTool::getDataWeight(const xAOD::EventInfo& eventInfo, const TString& trigger, bool mu_dependent) {
+
+ if(eventInfo.eventType(xAOD::EventInfo::IS_SIMULATION)) {
+ ATH_MSG_WARNING("Requesting Data Weight is not intended for simulated events. Returning 1.");
+ return 1; //no data weights for simulated events
+ }
+
+ if(!mu_dependent) return /*m_tool->*/m_activeTool->GetDataWeight(eventInfo.runNumber(), trigger);
+ return /*m_tool->*/m_activeTool->GetDataWeight( eventInfo.runNumber(), trigger, getLumiBlockMu(eventInfo) /*use the 'correct' mu instead of the one from the file!!*/ );
}
-} // namespace CP
\ No newline at end of file
+} // namespace CP
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/TPileupReweighting.cxx b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/TPileupReweighting.cxx
index cc6d00cc5f57f2da91b2d0b56c01e33d68a85e45..784d075a8b7cfc86b7eaade2225d3b249e244fc8 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/TPileupReweighting.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/Root/TPileupReweighting.cxx
@@ -36,189 +36,177 @@ Description: Tool to get the calculated MC pileup weight.
#include <TRandom3.h>
-ClassImp(Root::TPileupReweighting)
+ClassImp(CP::TPileupReweighting)
//=============================================================================
// Constructor
//=============================================================================
-Root::TPileupReweighting::TPileupReweighting(const char* name) :
+CP::TPileupReweighting::TPileupReweighting(const char* name) :
TNamed(name,"notitle"),
m_SetWarnings(true),m_debugging(false),
- m_countingMode(true),m_defaultChannel(0),m_unrepresentedDataAction(0),m_isInitialized(false),m_lumiVectorIsLoaded(false),
- m_dataScaleFactorX(1.),m_dataScaleFactorY(1.),m_dataScaleFactorZ(1.),
- m_mcScaleFactorX(1.),m_mcScaleFactorY(1.),m_mcScaleFactorZ(1.),
- m_nextPeriodNumber(1),m_ignoreFilePeriods(false),m_metadatatree(0),m_unrepDataTolerance(0.05),m_doGlobalDataWeight(false)
+ m_countingMode(true),m_unrepresentedDataAction(0),m_isInitialized(false),m_lumiVectorIsLoaded(false),
+ m_dataScaleFactorX(1.),m_dataScaleFactorY(1.),
+ m_mcScaleFactorX(1.),m_mcScaleFactorY(1.),
+ m_nextPeriodNumber(1),m_ignoreFilePeriods(false),m_metadatatree(0),m_unrepDataTolerance(0.05),m_doGlobalDataWeight(false),m_lumicalcRunNumberOffset(0), m_emptyHistogram(0), m_random3(0)
{
m_random3 = new TRandom3(0);
m_random3->SetSeed(1);
- //load the default pileup histogram binning
- //this is now done in the periodConfig (2012-05-09). We start with no default so that the default is picked up from the prw file
- //m_emptyHistograms["pileup"] = new TH1D("pileup_default","pileup_default",100,0,50);
- //m_emptyHistograms["pileup"]->SetDirectory(0);
- //set the lumi vector to size=11
- m_integratedLumiVector.resize(11,0.);
+
+ //create the global period
+ m_periods[-1] = new Period(-1,0,9999999,0/* the global default channel*/);
//MC12ab merger by default
- MergeMCRunNumbers(195847,195848);
+ //RemapPeriod(195848,195847);
+}
+
+
+
+//=============================================================================
+// Destructor
+//=============================================================================
+CP::TPileupReweighting::~TPileupReweighting() {
+ delete m_random3;
+}
+
+void CP::TPileupReweighting::RemapPeriod(Int_t periodNumber1, Int_t periodNumber2) {
+ //check if periodNumber2 exists
+ if(m_periods.find(periodNumber2)==m_periods.end()) {
+ m_periods[periodNumber2] = new Period(periodNumber2,0,0,m_periods[-1]->defaultChannel);
+ }
+ m_periods[periodNumber1] = m_periods[periodNumber2];
}
-Double_t Root::TPileupReweighting::GetIntegratedLumiFraction(Int_t mcRunNumber, UInt_t start, UInt_t end) {
- if(!m_isInitialized) {
- Error("GetIntegratedLumiFraction", "Please initialize the tool before retrieving lumi fractions");
- throw std::runtime_error("Throwing 1");
+void CP::TPileupReweighting::SetDefaultChannel(Int_t channel, Int_t periodNumber) {
+ //check if periodNumber2 exists
+ if(m_periods.find(periodNumber)==m_periods.end()) {
+ m_periods[periodNumber] = new Period(periodNumber,0,0,channel);
+ } else {
+ m_periods[periodNumber]->SetDefaultChannel(channel);
}
+}
+
+
+Int_t CP::TPileupReweighting::SetBinning(Int_t nbinsx, Double_t* xbins, Int_t nbinsy, Double_t* ybins) {
+ if(m_emptyHistogram) { delete m_emptyHistogram; }
+ if(nbinsy>0) m_emptyHistogram = new TH2D("default","default",nbinsx,xbins,nbinsy,ybins);
+ else m_emptyHistogram = new TH1D("default","default",nbinsx,xbins);
+ m_emptyHistogram->SetDirectory(0);
+ return 0;
+}
+Int_t CP::TPileupReweighting::SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy, Double_t ylow, Double_t yup) {
+ std::vector<double> xbins(nbinsx+1);
+ std::vector<double> ybins(nbinsy+1);
+ for(int i=0;i<nbinsx+1;i++) xbins[i] = xlow + i*(xup-xlow)/nbinsx;
+ for(int i=0;i<nbinsy+1;i++) ybins[i] = ylow + i*(yup-ylow)/nbinsy;
+ return SetBinning(nbinsx,&xbins[0],nbinsy,&ybins[0]);
+}
+Int_t CP::TPileupReweighting::SetBinning(TH1* hist) {
+ if(!hist) return -1;
+ if(m_emptyHistogram) { delete m_emptyHistogram; }
+ m_emptyHistogram = dynamic_cast<TH1*>(hist->Clone("default"));
+ m_emptyHistogram->SetDirectory(0);
+ return 0;
+}
+
+Int_t CP::TPileupReweighting::GetDefaultChannel(Int_t mcRunNumber) {
+ return m_periods[mcRunNumber]->defaultChannel;
+}
+
+
+Double_t CP::TPileupReweighting::GetIntegratedLumi(const TString& trigger) {
+ if(!m_isInitialized) { Info("GetIntegratedLumi","Initializing the subtool.."); Initialize(); }
if(!m_lumiVectorIsLoaded) {
- Error("GetIntegratedLumiFraction","No Lumicalc file loaded, so no lumi fraction possible, returning 0");
+ Error("GetIntegratedLumi","No Lumicalc file loaded, so cannot get integrated lumi, returning 0");
return 0;
}
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
- const std::map<Int_t, Double_t>::const_iterator it = periodTotals["pileup"][-1].find(mcRunNumber);
- if(it==periodTotals["pileup"][-1].end() || it->second == 0) {
- if(m_SetWarnings) Warning("GetIntegratedLumiFraction","%d has no associated lumi. Returning 0",mcRunNumber);
- return 0; /*throw 33;*/
+ if(trigger=="" || trigger=="None") return GetSumOfEventWeights(-1)/1E6;
+ //check the triggerPrimaryDistributions in the global period, if not present, build it
+ Period* global = m_periods[-1];
+ if(global->triggerHists.find(trigger)==global->triggerHists.end()) {
+ //try to do a subtrigger calculation we need to reopen all the lumicalc files and do the fiddily prescaled luminosity calculation
+ //will throw errors if can't
+ CalculatePrescaledLuminosityHistograms(trigger);
}
- //loop over assigned runs, if in range then include in numerator
- double total = 0;
- std::map<UInt_t, Double_t>::const_iterator itend = dataPeriodRunTotals["pileup"][mcRunNumber].end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = dataPeriodRunTotals["pileup"][mcRunNumber].begin(); it2 != itend; ++it2) {
- if(it2->first >= start && it2->first <= end) total += it2->second;
+ return global->triggerHists[trigger]->Integral(0,global->triggerHists[trigger]->GetNbinsX()+1)/1E6;
+}
+
+
+
+Double_t CP::TPileupReweighting::GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end) {
+ if(!m_isInitialized) { Info("GetIntegratedLumi","Initializing the subtool.."); Initialize(); }
+ //look through dataPeriodRunTotals["pileup"][periodNumber] for runs inside the given period
+ double total = 0;
+ for(auto run : m_periods[periodNumber]->runNumbers) {
+ if(run >= start && run <= end) total += m_runs[run].inputHists["None"]->GetSumOfWeights();
}
+ return total*1E-6;
+}
+
+Double_t CP::TPileupReweighting::GetLumiBlockIntegratedLumi(Int_t runNumber, UInt_t lb) {
+ if(!m_isInitialized) { Info("GetIntegratedLumi","Initializing the subtool.."); Initialize(); }
+ if(m_runs.find(runNumber)==m_runs.end()) return 0.;
+ auto& run = m_runs[runNumber];
- return total / it->second;
+ if(run.lumiByLbn.find(lb)==run.lumiByLbn.end()) return 0.;
+ return run.lumiByLbn[lb].first*1E-6;
}
-Double_t Root::TPileupReweighting::GetIntegratedLumiFraction(Int_t mcRunNumber, Double_t mu, UInt_t start, UInt_t end) {
- if(!m_isInitialized) {
- Error("GetIntegratedLumiFraction", "Please initialize the tool before retrieving lumi fractions");
- throw std::runtime_error("Throwing 1");
- }
+Float_t CP::TPileupReweighting::GetLumiBlockMu(Int_t runNumber, UInt_t lb) {
+ if(!m_isInitialized) { Info("GetLumiBlockMu","Initializing the subtool.."); Initialize(); }
+ if(m_runs.find(runNumber)==m_runs.end()) return 0.;
+ auto& run = m_runs[runNumber];
+
+ if(run.lumiByLbn.find(lb)==run.lumiByLbn.end()) return 0.;
+ return run.lumiByLbn[lb].second;
+}
+
+Double_t CP::TPileupReweighting::GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end) {
+ if(!m_isInitialized) { Info("GetIntegratedLumiFraction","Initializing the subtool.."); Initialize(); }
+
if(!m_lumiVectorIsLoaded) {
Error("GetIntegratedLumiFraction","No Lumicalc file loaded, so no lumi fraction possible, returning 0");
return 0;
}
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
-
- //determine the mubin
- TH1* hist = m_emptyHistograms["pileup"];
- if(!hist) { Error("GetIntegratedLumiFraction", "Cannot do this without a lumicalc file!"); return 0; }
- Int_t muBin = hist->FindFixBin(mu);
-
- //get the total lumi of this mcRunNumber with this mu - (was stored in the run=0 slot of the map)
- //check that this mcRunNumber actually has some lumi with this mu
- std::map<UInt_t, std::map<Int_t, std::map<UInt_t, Double_t> > >::iterator myMap1 = dataPeriodBinnedRunTotals["pileup"].find(muBin);
- if(myMap1 == dataPeriodBinnedRunTotals["pileup"].end()) {
- //no lumi with this mu
- return 0;
+ //total lumi in period is given by sumOfWeights[-1]
+ double total = m_periods[periodNumber]->sumOfWeights[-1];
+ //go through associated runs, adding up lumi from all within start and end (inclusive)
+ double numer(0);
+ for(auto run : m_periods[periodNumber]->runNumbers) {
+ if(run >= start && run <= end) numer += m_runs[run].inputHists["None"]->GetSumOfWeights();
}
- std::map<Int_t, std::map<UInt_t, Double_t> >::iterator myMap2 = myMap1->second.find(mcRunNumber);
- if(myMap2 == myMap1->second.end()) {
- //no lumi with this mu in this mcRunNumber
+ return numer / total;
+}
+
+Double_t CP::TPileupReweighting::GetIntegratedLumiFraction(Int_t periodNumber, Double_t mu, UInt_t start, UInt_t end) {
+ if(!m_isInitialized) { Info("GetIntegratedLumiFraction","Initializing the subtool.."); Initialize(); }
+ if(!m_lumiVectorIsLoaded) {
+ Error("GetIntegratedLumiFraction","No Lumicalc file loaded, so no lumi fraction possible, returning 0");
return 0;
}
+ //determine the mubin
+ if(!m_emptyHistogram) { Error("GetIntegratedLumiFraction", "Cannot do this without a lumicalc file!"); return 0; }
+ Int_t muBin = m_emptyHistogram->FindFixBin(mu);
+ //the total mu in this bin in this period is given by the triggerHist = "None" hist
+ double total = m_periods[periodNumber]->triggerHists["None"]->GetBinContent(muBin);
//loop over assigned runs, if in range then include in numerator
- double total = 0;
- std::map<UInt_t, Double_t>::const_iterator itend = myMap2->second.end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = myMap2->second.begin(); it2 != itend; ++it2) {
- if(it2->first >= start && it2->first <= end) total += it2->second;
+ double numer = 0;
+ for(auto run : m_periods[periodNumber]->runNumbers) {
+ if(run >= start && run <= end) numer += m_runs[run].muDist->GetBinContent(muBin);
}
- return total / myMap2->second[0];
-}
-
-Double_t Root::TPileupReweighting::GetIntegratedLumi(UInt_t start, UInt_t end) {
- //look through dataPeriodRunTotals["pileup"][-1] for runs inside the given period
- double total = 0;
- std::map<UInt_t, Double_t>::const_iterator itend = dataPeriodRunTotals["pileup"][-1].end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = dataPeriodRunTotals["pileup"][-1].begin(); it2 != itend; ++it2) {
- if(it2->first >= start && it2->first <= end) total += it2->second;
- }
- return total/1E6;
+ return numer / total;
}
-Double_t Root::TPileupReweighting::GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end) {
- //look through dataPeriodRunTotals["pileup"][periodNumber] for runs inside the given period
- double total = 0;
- std::map<UInt_t, Double_t>::const_iterator itend = dataPeriodRunTotals["pileup"][periodNumber].end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = dataPeriodRunTotals["pileup"][periodNumber].begin(); it2 != itend; ++it2) {
- if(it2->first >= start && it2->first <= end) total += it2->second;
- }
- return total/1E6;
-}
-//=============================================================================
-// Destructor
-//=============================================================================
-Root::TPileupReweighting::~TPileupReweighting() {
- delete m_random3;
- //delete all histograms
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > >::iterator it0=m_inputHistograms.begin();it0!=m_inputHistograms.end();++it0) {
- for(std::map<Int_t,std::map<Int_t, TH1*> >::iterator it1=it0->second.begin();it1!=it0->second.end();++it1) {
- for(std::map<Int_t, TH1*>::iterator it2=it1->second.begin();it2!=it1->second.end();++it2) {
- if(it2->second) delete it2->second;
- }
- it1->second.clear();
- }
- it0->second.clear();
- }
- m_inputHistograms.clear();
- if(m_debugging) Info("destructor","Deleting primary distributions");
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH1D*> > >::iterator it0=primaryDistributions.begin();it0!=primaryDistributions.end();++it0) {
- for(std::map<Int_t,std::map<Int_t, TH1D*> >::iterator it1=it0->second.begin();it1!=it0->second.end();++it1) {
- for(std::map<Int_t, TH1D*>::iterator it2=it1->second.begin();it2!=it1->second.end();++it2) {
- if(it2->second) {
- if(m_debugging) { Info("destructor","weight=%s chan=%d run=%d",it0->first.Data(),it1->first,it2->first);
- Info("destructor","Deleting %s",it2->second->GetName());
- }
- delete it2->second;
- }
- }
- it1->second.clear();
- }
- it0->second.clear();
- }
- primaryDistributions.clear();
- if(m_debugging) Info("destructor","Deleting secondary distributions");
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH2D*> > >::iterator it0=secondaryDistributions.begin();it0!=secondaryDistributions.end();++it0) {
- for(std::map<Int_t,std::map<Int_t, TH2D*> >::iterator it1=it0->second.begin();it1!=it0->second.end();++it1) {
- for(std::map<Int_t, TH2D*>::iterator it2=it1->second.begin();it2!=it1->second.end();++it2) {
- if(it2->second) {
- if(m_debugging) { Info("destructor","weight=%s chan=%d run=%d",it0->first.Data(),it1->first,it2->first);
- Info("destructor","Deleting %s",it2->second->GetName());
- }
- delete it2->second;
- }
- }
- it1->second.clear();
- }
- it0->second.clear();
- }
- secondaryDistributions.clear();
- if(m_debugging) Info("destructor","Deleting trigger primary distributions");
- for(std::map<TString, std::map<Int_t, TH1D*> >::iterator it0=m_triggerPrimaryDistributions.begin();it0!=m_triggerPrimaryDistributions.end();++it0) {
- for(std::map<Int_t, TH1D*>::iterator it=it0->second.begin();it!=it0->second.end();++it) {
- if(it->second) {
- if(m_debugging) { Info("destructor","trigger=%s run=%d",it0->first.Data(),it->first);
- Info("destructor","Deleting %s",it->second->GetName());
- }
- delete it->second;
- }
- }
- it0->second.clear();
- }
- m_triggerPrimaryDistributions.clear();
-
-}
-Int_t Root::TPileupReweighting::UsePeriodConfig(const TString configName) {
+Int_t CP::TPileupReweighting::UsePeriodConfig(const TString configName) {
if(configName=="MC11a") {
AddPeriod(180164, 177986,180481); //associates mc runnumber 180164 with data period 177986 to 180481 (period B-D)
AddPeriod(183003, 180614,184169); //period E-H
AddPeriod(185649, 185353,186934); //period I-K1. For I-K you would change the last number to 187815
AddPeriod(185761, 186935,191933); //everything else. Thanks Ellie!
- AddBinning("pileup",100,0,50);
+ SetUniformBinning(100,0,50);
Info("UsePeriodConfig","Using MC11a Period configuration");
return 0;
} else if(configName=="MC11b" || configName=="MC11c") {
@@ -226,192 +214,150 @@ Int_t Root::TPileupReweighting::UsePeriodConfig(const TString configName) {
AddPeriod(183003, 180614, 184169);
AddPeriod(186169, 185353, 187815);
AddPeriod(189751, 188902, 191933);
- AddBinning("pileup",100,0,50);
+ SetUniformBinning(100,0,50);
Info("UsePeriodConfig","Using MC11b/c Period configuration");
return 0;
} else if(configName=="MC12a") {
AddPeriod(195847,200804,216432);
//mc12a binning is in integer values of mu
- AddBinning("pileup",50,-0.5,49.5);
+
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC12a, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
+ return -1;
+ }
+
+ SetUniformBinning(50,-0.5,49.5);
Info("UsePeriodConfig","Using MC12a Period configuration");
return 0;
} else if(configName=="MC12b") {
AddPeriod(195848,200804,216432);
//mc12b binning is in integer values of mu
- AddBinning("pileup",50,-0.5,49.5);
+
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC12b, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
+ return -1;
+ }
+
+ SetUniformBinning(50,-0.5,49.5);
Info("UsePeriodConfig","Using MC12b Period configuration");
return 0;
} else if(configName=="MC12ab") {
AddPeriod(195847,200804,216432);
- //disable warnings temporarily for next AddPeriod, as we know we are causing a conflict
- bool tmp = m_SetWarnings;
- m_SetWarnings=false;
AddPeriod(195848,200804,216432);
- m_SetWarnings = tmp;
+
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC12ab, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
+ return -1;
+ }
+
//mc12a/b binning is in integer values of mu
- AddBinning("pileup",50,-0.5,49.5);
+ SetUniformBinning(50,-0.5,49.5);
Info("UsePeriodConfig","Using MC12ab Period configuration");
return 0;
- } else if(configName=="MC14") {
+ } else if(configName=="MC14_8TeV") {
AddPeriod(212272,200804,216432);
- AddBinning("pileup",50,-0.5,49.5);
- Info("UsePeriodConfig","Using MC14 Period configuration");
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01 || fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC14_8TeV, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC14_8TeV, an incompatible config has already been set up");
+ return -1;
+ }
+ SetUniformBinning(50,-0.5,49.5);
+ Info("UsePeriodConfig","Using MC14_8TeV Period configuration");
+ return 0;
+ } else if(configName=="MC14_13TeV") {
+ AddPeriod(222222,222222,999999);
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=100 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(100)-99.5)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC14_13TeV, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC14_13TeV, an incompatible config has already been set up");
+ return -1;
+ }
+ SetUniformBinning(100,-0.5,99.5);
+ Info("UsePeriodConfig","Using MC14_13TeV Period configuration");
+ return 0;
+ } else if(configName=="MC15") {
+ AddPeriod(222510,222222,999999);
+ AddPeriod(222525,222222,999999);
+ if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=100 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1))>0.01 || fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(100)-100.)>0.01) ) {
+ Error("UsePeriodConfig","Cannot use MC15, an incompatible config has already been set up");
+ throw std::runtime_error("Throwing 13: Cannot use MC15, an incompatible config has already been set up");
+ return -1;
+ }
+ SetUniformBinning(100,0,100); //Thanks Eric </sarcasm>
+ Info("UsePeriodConfig","Using MC15 Period configuration");
return 0;
}
Error("UsePeriodConfig","Unrecognized period config");
return -1;
}
-std::vector<Double_t> Root::TPileupReweighting::getIntegratedLumiVector() {
- if(!m_lumiVectorIsLoaded) {
- Error("getIntegratedLumiVector","The vector has not been filled. You must AddDataDistribution, using the iLumiCalc file, to fill it");
- throw std::runtime_error("Throwing 99");
- }
- if(!m_isInitialized) {
- Error("getIntegratedLumiVector","Please initialize the tool before retrieving the lumi vector");
- throw std::runtime_error("Throwing 1");
- }
- return m_integratedLumiVector;
-}
///returns a PeriodID. These count up from 1
-Int_t Root::TPileupReweighting::AddPeriod(Int_t mcRunNumber, UInt_t start, UInt_t end) {
- if(m_isInitialized) {
- Error("AddPeriod","You cannot AddPeriod after initializing the tool. Reorder your code!");
- throw std::runtime_error("Throwing 1: You cannot AddPeriod after initializing the tool. Reorder your code!");
- }
- //check if this period is already defined
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=m_periods[mcRunNumber].begin();periods!=m_periods[mcRunNumber].end();++periods) {
- std::pair<UInt_t,UInt_t>& period = periods->second;
- if(period.first == start && period.second==end) {
- return 0;
- } else if(start>=period.first && end<=period.second) {
- //fully contained in existing period, just print a warning
- Warning("AddPeriod","Period (%d,%d,%d) fully contained within existing period (%d,%d,%d). Doing nothing.",mcRunNumber,start,end,mcRunNumber,period.first,period.second);
- return 0;
- } else if(start==period.first && end==999999) {
- //special case where we are 'extending' the period to infinity. This is a kludge to make old config files hadded to new ones work
- Warning("AddPeriod","Extending period [%d,%d] to infinity. I hope this only happened because you are using an old prw config file (made pre 00-02-10)",period.first,period.second);
- period.second = 999999;
- return 0;
- } else if((start <= period.second && start >= period.first) || (end <= period.second && end >= period.first) ) {
- //overlaps with existing period. reject this Add
- Error("AddPeriod","(%d,%d,%d) overlaps with existing period setup",mcRunNumber,start,end);
- throw std::runtime_error("Throwing 2: (%d,%d,%d) overlaps with existing period setup");
- }
- }
- //check if the period being added conflicts with a period in another mcRunNumber... we just print a warning at this stage... not sure if it is a problem...
- for(std::map<Int_t, std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator periodRuns=m_periods.begin();periodRuns!=m_periods.end();++periodRuns) {
- if(periodRuns->first==mcRunNumber) continue; //already checked this one
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=periodRuns->second.begin();periods!=periodRuns->second.end();++periods) {
- std::pair<UInt_t,UInt_t>& period = periods->second;
- //if the period start and end match perfectly and this is a registered mc merger, then no warning
- if(start==period.first && end==period.second && m_mcMerges.find(periodRuns->first)!=m_mcMerges.end()&& m_mcMerges[periodRuns->first].find(periodRuns->first)!=m_mcMerges[periodRuns->first].end() && m_mcMerges[periodRuns->first][periodRuns->first]==mcRunNumber) continue;
- if(start==period.first && end==period.second && m_mcMerges.find(mcRunNumber)!=m_mcMerges.end()&& m_mcMerges[mcRunNumber].find(mcRunNumber)!=m_mcMerges[mcRunNumber].end() && m_mcMerges[mcRunNumber][mcRunNumber]==periodRuns->first) continue;
- if((start <= period.second && start >= period.first) || (end <= period.second && end >= period.first) ) {
- //overlaps with existing period. show a warning
- if(m_SetWarnings) Warning("AddPeriod","(%d,%d,%d) conflicts with existing period setup (%d,%d,%d)",mcRunNumber,start,end,periodRuns->first,period.first,period.second);
- }
- }
- }
- std::pair<UInt_t,UInt_t> thePair; thePair.first=start;thePair.second=end;
- m_periods[mcRunNumber][m_nextPeriodNumber] = thePair;
- //also put this in a handy backmapping for getting this period's periodWeight and mc distributions (uses same ones)
- m_periodToMCRun[m_nextPeriodNumber]=mcRunNumber;
- m_nextPeriodNumber++;
- return m_nextPeriodNumber-1;
-}
+Int_t CP::TPileupReweighting::AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end) {
+ //you can add more periods only in counting mode, but they must not be with subperiods!
-Int_t Root::TPileupReweighting::GetPeriodNumber(Int_t mcRunNumber, UInt_t start, UInt_t end) {
- //try to locate this period number in the map
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=m_periods[mcRunNumber].begin();periods!=m_periods[mcRunNumber].end();++periods) {
- std::pair<UInt_t,UInt_t> period = periods->second;
- if(period.first == start && period.second==end) {
- return periods->first;
- }
- }
- return -1;
-}
+ if(m_isInitialized && !m_countingMode) {
+ Error("AddPeriod","You cannot AddPeriod after initializing the tool, except when in config file generating mode. Reorder your code!");
+ throw std::runtime_error("Throwing 1: You cannot AddPeriod after initializing the tool, except when in config file generating mode. Reorder your code!");
+ }
-Int_t Root::TPileupReweighting::GetFirstFoundPeriodNumber(UInt_t runNumber) {
- //try to locate this period number in the map
- for(std::map<Int_t,std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator it = m_periods.begin();it!=m_periods.end();++it) {
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=it->second.begin();periods!=it->second.end();++periods) {
- std::pair<UInt_t,UInt_t> period = periods->second;
- if(period.first <= runNumber && period.second>=runNumber) {
- return periods->first;
- }
- }
+ if(m_periods.find(periodNumber)==m_periods.end()) {
+ m_periods[periodNumber] = new Period(periodNumber,start,end,m_periods[-1]->defaultChannel);
+ return periodNumber;
}
- return -1;
-}
-Int_t Root::TPileupReweighting::GetPeriodNumber(Int_t mcRunNumber, UInt_t runNumber) {
- //try to locate this period number in the map
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=m_periods[mcRunNumber].begin();periods!=m_periods[mcRunNumber].end();++periods) {
- std::pair<UInt_t,UInt_t> period = periods->second;
- if(period.first <= runNumber && period.second>=runNumber) {
- return periods->first;
- }
- }
- return -1;
-}
+ //if period exists but both start and end are zero, then just assign that
+ Period* p = m_periods[periodNumber];
-Int_t Root::TPileupReweighting::AddBinning(const TString weightName,Int_t nbinsx, Double_t* xbins) {
- if(m_emptyHistograms.find(weightName)!=m_emptyHistograms.end()) {
- delete m_emptyHistograms[weightName];
- }
- m_emptyHistograms[weightName] = new TH1D(weightName+"_default","default",nbinsx,xbins);
- m_emptyHistograms[weightName]->SetDirectory(0);
- return 0;
-}
-Int_t Root::TPileupReweighting::AddBinning(const TString weightName,Int_t nbinsx, Double_t xlow, Double_t xup) {
- if(m_emptyHistograms.find(weightName)!=m_emptyHistograms.end()) {
- delete m_emptyHistograms[weightName];
+ if(p->subPeriods.size() == 0 && p->start==0 && p->end==0) { p->start=start; p->end=end; return periodNumber; }
+
+ //check if period already contained (i.e. same period or one of the existing subperiods) ... if it is, do nothing
+ if(p->start==start && p->end==end) return periodNumber;
+ for(auto pp : p->subPeriods) {
+ if(pp->start==start && pp->end==end) return pp->id;
}
- m_emptyHistograms[weightName] = new TH1D(weightName+"_default","default",nbinsx,xlow,xup);
- m_emptyHistograms[weightName]->SetDirectory(0);
- return 0;
-}
-Int_t Root::TPileupReweighting::AddBinning(const TString weightName,Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy, Double_t ylow, Double_t yup) {
- if(m_emptyHistograms.find(weightName)!=m_emptyHistograms.end()) {
- delete m_emptyHistograms[weightName];
+
+ //if get here, and in counting mode, this isn't allowed! .. i.e. we can't have subperiods
+ if(m_isInitialized && m_countingMode) {
+ Error("AddPeriod","You cannot have subperiods when in Config File Generating mode");
+ throw std::runtime_error("Throwing 44: In Config File Generating mode, but detected subperiods in the period definition. This is not allowed.");
}
- m_emptyHistograms[weightName] = new TH2D(weightName+"_default","default",nbinsx,xlow,xup,nbinsy,ylow,yup);
- m_emptyHistograms[weightName]->SetDirectory(0);
- return 0;
-}
-Int_t Root::TPileupReweighting::AddBinning(const TString weightName,TH1* hist) {
- if(!hist) return -1;
- if(m_emptyHistograms.find(weightName)!=m_emptyHistograms.end()) {
- delete m_emptyHistograms[weightName];
+ //if period is already defined, we create a new period with start and end, and assign our period to this one
+ //also create a new period for the existing period configuration of this period
+
+ if(p->subPeriods.size()==0) {
+ while(m_periods.find(m_nextPeriodNumber) != m_periods.end()) m_nextPeriodNumber++;
+ m_periods[m_nextPeriodNumber] = new Period(m_nextPeriodNumber,p->start,p->end,p->defaultChannel);
+ p->subPeriods.push_back(m_periods[m_nextPeriodNumber]);
+ p->start = 0; p->end=0;
}
- m_emptyHistograms[weightName] = dynamic_cast<TH1*>(hist->Clone(weightName+"_default"));
- m_emptyHistograms[weightName]->SetDirectory(0);
- return 0;
-}
-Int_t Root::TPileupReweighting::GetDefaultChannel(Int_t mcRunNumber) {
- if(m_defaultChannelsByMCRunNumber.find(mcRunNumber)!=m_defaultChannelsByMCRunNumber.end()) return m_defaultChannelsByMCRunNumber[mcRunNumber];
- else return m_defaultChannel;
+ while(m_periods.find(m_nextPeriodNumber) != m_periods.end()) m_nextPeriodNumber++;
+ m_periods[m_nextPeriodNumber] = new Period(m_nextPeriodNumber,start,end, p->defaultChannel);
+ p->subPeriods.push_back(m_periods[m_nextPeriodNumber]);
+
+ return m_nextPeriodNumber;
}
-//combine run number's histograms, periods, everything, then remove trace of second run number
-void Root::TPileupReweighting::MergeMCRunNumbers(Int_t mcRunNumber1, Int_t mcRunNumber2, Int_t newMCRunNumber) {
- if(m_isInitialized) {
- Error("MergeMCRunNumbers","You cannot MergeMCRunNumbers after initializing the tool. Reorder your code!");
- throw std::runtime_error("Throwing 5: You cannot MergeMCRunNumbers after initializing the tool. Reorder your code!");
+Int_t CP::TPileupReweighting::GetFirstFoundPeriodNumber(UInt_t runNumber) {
+ //go through periods, first period that includes this runNumber wins!
+ //otherwise it returns the 'global period number': -1
+ for(auto period : m_periods) {
+ if(period.second->id==-1) continue;
+ if(period.second->contains(runNumber)) return period.second->id;
}
- if(newMCRunNumber==0) newMCRunNumber=mcRunNumber1;
- m_mcMerges[newMCRunNumber][mcRunNumber1]=mcRunNumber2;
+ return -1;
}
-TH1* Root::TPileupReweighting::CloneEmptyHistogram(const TString weightName, Int_t runNumber, Int_t channelNumber) {
- TString s(weightName);
+TH1* CP::TPileupReweighting::CloneEmptyHistogram(Int_t runNumber, Int_t channelNumber) {
+
+ TString s("pileup");
if(channelNumber>=0) {
s += "_chan"; s += channelNumber;
} else {
@@ -419,12 +365,12 @@ TH1* Root::TPileupReweighting::CloneEmptyHistogram(const TString weightName, Int
}
s+="_run"; s+= runNumber;
//get the empty histogram
- if(m_emptyHistograms.find(weightName)==m_emptyHistograms.end()) {
- Error("CloneEmptyHistogram","There is no binning info for weight %s - use AddBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)",weightName.Data());
- throw std::runtime_error("Throwing 47: There is no binning info - use AddBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)");
+ if(!m_emptyHistogram) {
+ Error("CloneEmptyHistogram","There is no binning info - use SetBinning/SetUniformBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)");
+ throw std::runtime_error("Throwing 47: There is no binning info - use SetBinning/SetUniformBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)");
}
- TH1* out = dynamic_cast<TH1*>(m_emptyHistograms[weightName]->Clone(s));
+ TH1* out = dynamic_cast<TH1*>(m_emptyHistogram->Clone(s));
out->SetTitle(s);
out->SetDirectory(0); //take ownership
out->SetEntries(0);
@@ -432,7 +378,7 @@ TH1* Root::TPileupReweighting::CloneEmptyHistogram(const TString weightName, Int
}
-Int_t Root::TPileupReweighting::GenerateMetaDataFile(const TString fileName,const TString channelBranchName) {
+Int_t CP::TPileupReweighting::GenerateMetaDataFile(const TString fileName,const TString channelBranchName) {
TTree inTree("in","in");
inTree.ReadFile(fileName);
@@ -476,7 +422,7 @@ Int_t Root::TPileupReweighting::GenerateMetaDataFile(const TString fileName,cons
return 0;
}
-Int_t Root::TPileupReweighting::AddMetaDataFile(const TString fileName,const TString channelBranchName) {
+Int_t CP::TPileupReweighting::AddMetaDataFile(const TString fileName,const TString channelBranchName) {
TDirectory* origDir = gDirectory;
TTree* tmp = 0;
TFile* rootFile = 0;
@@ -533,7 +479,7 @@ Int_t Root::TPileupReweighting::AddMetaDataFile(const TString fileName,const TSt
}
}
- if(tmp) delete tmp;
+ delete tmp;
if(rootFile) { rootFile->Close();delete rootFile;}
// Return to the original ROOT directory
@@ -542,7 +488,7 @@ Int_t Root::TPileupReweighting::AddMetaDataFile(const TString fileName,const TSt
return 0;
}
-TTree* Root::TPileupReweighting::GetMetaDataTree() {
+TTree* CP::TPileupReweighting::GetMetaDataTree() {
if(m_metadatatree) delete m_metadatatree;
m_metadatatree = new TTree(TString(this->GetName())+"MetaData",TString(this->GetName())+"MetaData");
m_metadatatree->SetDirectory(0);
@@ -567,7 +513,8 @@ TTree* Root::TPileupReweighting::GetMetaDataTree() {
data["SumOfEventWeights"]=0.;
m_metadatatree->Branch("SumOfEventWeights",&(data["SumOfEventWeights"]));
//and add to channels list any channels that have events
- for(std::map<Int_t,Double_t>::iterator it=globalNumberOfEntries["pileup"].begin();it!=globalNumberOfEntries["pileup"].end();++it) {
+ auto& myMap = m_periods[-1]->sumOfWeights;
+ for(std::map<Int_t,Double_t>::iterator it=myMap.begin();it!=myMap.end();++it) {
if(it->first>=0 && it->second>0.) channels[it->first]=true;
}
@@ -577,27 +524,29 @@ TTree* Root::TPileupReweighting::GetMetaDataTree() {
for(std::map<TString,Double_t>::iterator it2=data.begin();it2!=data.end();++it2) {
if(it2->first=="NumberOfEvents") {
//check for this in globalNumberOfEntries
- if(globalNumberOfEntries["pileup"].find(channel)==globalNumberOfEntries["pileup"].end()) {
+ if(myMap.find(channel)==myMap.end()) {
data[it2->first]=0.;
Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
} else {
- data[it2->first]=globalNumberOfEntries["pileup"][channel];
+ data[it2->first]=m_periods[-1]->numberOfEntries[channel];
}
} else if(it2->first=="SumOfEventWeights") {
//check for this in globalTotals
- if(globalTotals["pileup"].find(channel)==globalTotals["pileup"].end()) {
+ auto& myMap2 = m_periods[-1]->sumOfWeights;
+ if(myMap2.find(channel)==myMap2.end()) {
data[it2->first]=0.;
Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
} else {
- data[it2->first]=globalTotals["pileup"][channel];
+ data[it2->first]=m_periods[-1]->sumOfWeights[channel];
}
} else {
- if(m_metadata[it2->first].find(channel)==m_metadata[it2->first].end()) {
+ auto& myMap2 = m_metadata[it2->first];
+ if(myMap2.find(channel)==myMap2.end()) {
//this channel doesn't have this property
data[it2->first]=0.;
Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
} else {
- data[it2->first]=m_metadata[it2->first][channel];
+ data[it2->first]=myMap2[channel];
}
}
}
@@ -610,11 +559,14 @@ TTree* Root::TPileupReweighting::GetMetaDataTree() {
return m_metadatatree;
}
-Int_t Root::TPileupReweighting::AddDistribution(TH1* hist, Int_t runNumber, Int_t channelNumber) {
- return AddDistribution(hist,"pileup",runNumber,channelNumber);
-}
-void Root::TPileupReweighting::AddDistributionTree(TTree *tree, TFile *file) {
+
+
+
+
+
+
+void CP::TPileupReweighting::AddDistributionTree(TTree *tree, TFile *file) {
bool isMC=true;
//should have the standard structure
@@ -660,9 +612,6 @@ void Root::TPileupReweighting::AddDistributionTree(TTree *tree, TFile *file) {
tree->GetEntry(i);
TString sHistName(histName);
TString weightName(customName);
- if(hasDefaultsBranch) {
- if(isDefaultForRunNumber) m_defaultChannelsByMCRunNumber[runNbr]=channel;
- }
if(loadedHistos.find(sHistName)==loadedHistos.end()) {
loadedHistos[sHistName]=true;
if(!m_ignoreFilePeriods && isMC) {
@@ -673,145 +622,97 @@ void Root::TPileupReweighting::AddDistributionTree(TTree *tree, TFile *file) {
}
//load the histogram
TH1 *histo = (TH1*)file->Get( sHistName );
+ if(!histo) histo = (TH1*)file->Get( TString::Format("%sPileupReweighting/%s",m_prwFilesPathPrefix.c_str(),sHistName.Data()) ); //try from subdir
if(!histo) {
Error("AddDistributionTree","Unable to find the histogram %s in the File %s",sHistName.Data(),file->GetName());
throw std::runtime_error("Throwing 21");
}
//add it to the histograms
- AddDistribution(histo, weightName, runNbr, channel);
+ AddDistribution(histo, runNbr, channel);
+ //we also add it to the global count, if this isn't data
+ if(channel>=0) AddDistribution(histo,-1,channel);
+ }
+ if(hasDefaultsBranch) {
+ if(isDefaultForRunNumber && m_periods.find(runNbr)!=m_periods.end()) m_periods[runNbr]->defaultChannel=channel;
}
}
}
//data is signalled by a negative channel number
-Int_t Root::TPileupReweighting::AddDistribution(TH1* hist,const TString weightName, Int_t runNumber, Int_t channelNumber) {
+Int_t CP::TPileupReweighting::AddDistribution(TH1* hist,Int_t runNumber, Int_t channelNumber) {
if(m_isInitialized) {
Error("AddDistribution","You cannot AddDistribution after initializing the tool. Reorder your code!");
throw std::runtime_error("Throwing 5: You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
}
-
- //std::map<Int_t,TH1*>& inputHistograms = (channelNumber>=0) ? m_inputMCHistograms[weightName][channelNumber] : m_inputDataHistograms[weightName];
- std::map<Int_t,TH1*>& inputHistograms = m_inputHistograms[weightName][channelNumber];
+ if(channelNumber>=0 && !m_periods[runNumber]) {
+ Error("AddDistribution","Unrecognised periodNumber: %d .. please use AddPeriod to define a period",runNumber);
+ throw std::runtime_error("Throwing 6: Unrecognised periodNumber. Please use AddPeriod to define a period");
+ }
- //check if we need to create an empty histogram
- if(inputHistograms.find(runNumber) == inputHistograms.end()) {
+ TH1*& inputHist = (channelNumber<0) ? m_runs[runNumber].inputHists["None"] : m_periods[runNumber]->inputHists[channelNumber];
+
+ if(!inputHist) {
//if no emptyHistogram specified, we will use this histogram as the structure for the empty;
- if(m_emptyHistograms.find(weightName)==m_emptyHistograms.end()) {
- TString s = weightName; s+= "_default";
+ if(!m_emptyHistogram) {
//check if the input histogram is a TH1D or not. We like TH1D, not yucky TH1F
if(strcmp(hist->IsA()->GetName(),"TH1D")==0 || strcmp(hist->IsA()->GetName(),"TH2D")==0) {
- m_emptyHistograms[weightName] = dynamic_cast<TH1*>(hist->Clone(s));
+ m_emptyHistogram = dynamic_cast<TH1*>(hist->Clone("default"));
} else {
//check dimensionality
if(hist->GetDimension()==1) {
std::vector<Double_t> binsX;
for(int i=0;i<=hist->GetNbinsX();i++) binsX.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
TH1D tmpHist("tmpHist","tmpHist",binsX.size()-1,&binsX[0]);
- m_emptyHistograms[weightName] = dynamic_cast<TH1*>(tmpHist.Clone(s));
+ m_emptyHistogram = dynamic_cast<TH1*>(tmpHist.Clone("default"));
} else if(hist->GetDimension()==2) {
std::vector<Double_t> binsX;std::vector<Double_t> binsY;
for(int i=0;i<=hist->GetNbinsX();i++) binsX.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
for(int i=0;i<=hist->GetNbinsY();i++) binsY.push_back(hist->GetYaxis()->GetBinLowEdge(i+1));
TH2D tmpHist("tmpHist","tmpHist",binsX.size()-1,&binsX[0],binsY.size()-1,&binsY[0]);
- m_emptyHistograms[weightName] = dynamic_cast<TH1*>(tmpHist.Clone(s));
- } else if(hist->GetDimension()==3) {
- std::vector<Double_t> binsX;std::vector<Double_t> binsY;std::vector<Double_t> binsZ;
- for(int i=0;i<=hist->GetNbinsX();i++) binsX.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
- for(int i=0;i<=hist->GetNbinsY();i++) binsY.push_back(hist->GetYaxis()->GetBinLowEdge(i+1));
- for(int i=0;i<=hist->GetNbinsZ();i++) binsZ.push_back(hist->GetZaxis()->GetBinLowEdge(i+1));
- TH3D tmpHist("tmpHist","tmpHist",binsX.size()-1,&binsX[0],binsY.size()-1,&binsY[0],binsZ.size()-1,&binsZ[0]);
- m_emptyHistograms[weightName] = dynamic_cast<TH1*>(tmpHist.Clone(s));
+ m_emptyHistogram = dynamic_cast<TH1*>(tmpHist.Clone("default"));
} else {
Error("AddDistribution","Unknown input histogram dimensionality: %d",hist->GetDimension());
throw std::runtime_error("Throwing 98");
}
}
- m_emptyHistograms[weightName]->Reset();//clear the histogram
- m_emptyHistograms[weightName]->SetEntries(0);
- m_emptyHistograms[weightName]->SetDirectory(0);
+ m_emptyHistogram->Reset();m_emptyHistogram->SetEntries(0);m_emptyHistogram->SetDirectory(0);//clear the histogram
}
- inputHistograms[runNumber] = CloneEmptyHistogram(weightName,runNumber,channelNumber);
+ inputHist = CloneEmptyHistogram(runNumber,channelNumber);
+
}
+
+
//iterator over bins of histogram, filling the TH1 stored in the data map
- Double_t numEntries = inputHistograms[runNumber]->GetEntries();
- Int_t bin,binx,biny,binz;
- for(binz=1; binz<=hist->GetNbinsZ(); binz++) {
- for(biny=1; biny<=hist->GetNbinsY(); biny++) {
- for(binx=1; binx<=hist->GetNbinsX(); binx++) {
- bin = hist->GetBin(binx,biny,binz);
- Double_t value = hist->GetBinContent(bin);
- Double_t x = hist->GetXaxis()->GetBinCenter(binx);
- Double_t y = hist->GetYaxis()->GetBinCenter(biny);
- Double_t z = hist->GetZaxis()->GetBinCenter(binz);
- //shift x,y,z by the MCScaleFactors as appropriate
- if(channelNumber>=0) {x *= m_mcScaleFactorX; y *= m_mcScaleFactorY; z *= m_mcScaleFactorZ;}
- else { x *= m_dataScaleFactorX; y *= m_dataScaleFactorY; z *= m_dataScaleFactorZ; }
-#if ROOT_VERSION_CODE >= ROOT_VERSION(5,28,0)
- Int_t inBin = inputHistograms[runNumber]->FindFixBin(x,y,z);
-#else
- Int_t inBin = inputHistograms[runNumber]->GetXaxis()->FindFixBin(x);
-#endif
- Double_t inValue = inputHistograms[runNumber]->GetBinContent(inBin);
- inputHistograms[runNumber]->SetBinContent(inBin,inValue+value);
- }
+ Double_t numEntries = inputHist->GetEntries();
+ Int_t bin,binx,biny;
+ for(biny=1; biny<=hist->GetNbinsY(); biny++) {
+ for(binx=1; binx<=hist->GetNbinsX(); binx++) {
+ bin = hist->GetBin(binx,biny);
+ Double_t value = hist->GetBinContent(bin);
+ Double_t x = hist->GetXaxis()->GetBinCenter(binx);
+ Double_t y = hist->GetYaxis()->GetBinCenter(biny);
+ //shift x,y,z by the MCScaleFactors as appropriate
+ if(channelNumber>=0) {x *= m_mcScaleFactorX; y *= m_mcScaleFactorY;}
+ else { x *= m_dataScaleFactorX; y *= m_dataScaleFactorY; }
+ Int_t inBin = inputHist->FindFixBin(x,y);
+ Double_t inValue = inputHist->GetBinContent(inBin);
+ inputHist->SetBinContent(inBin,inValue+value);
}
}
+
//also keep track of the number of entries
//SetBinContent screws with the entry count, so had to record it before the loops above
- inputHistograms[runNumber]->SetEntries(numEntries+hist->GetEntries());
+ inputHist->SetEntries(numEntries+hist->GetEntries());
m_countingMode=false;
return 0;
}
-Int_t Root::TPileupReweighting::AddDistributionFromFile(const TString fileName, const TString histName, const TString weightName, Int_t runNumber, Int_t channelNumber) {
- if(m_isInitialized) {
- Error("AddDistributionFromFile","You cannot AddDistributionFromFile after initializing the tool. Reorder your code!");
- throw std::runtime_error("Throwing 5: You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
- }
-
- // Cache the current directory in the root file
- TDirectory* origDir = gDirectory;
- // Load the data ROOT file
- TFile* rootFile = TFile::Open( fileName, "READ" );
- if ( rootFile->IsZombie() ) {
- Error("AddConfigFile","Could not open file: %s",fileName.Data());
- std::string toThrow = "Throwing 6: Could not open file: "; toThrow += fileName;
- throw std::runtime_error(toThrow);
- } else {
- //decide if this is a histogram or a ttree
- if(!rootFile->Get(histName)) {
- Error("AddDistributionFromFile","Could not find object %s in file %s",histName.Data(),fileName.Data());
- throw std::runtime_error("Throwing 7");
- }
- if(rootFile->Get(histName)->InheritsFrom("TH1")) {
- TH1 *tmp = (TH1*)rootFile->Get( histName );
- //tmp->SetDirectory(0);//take ownership
- AddDistribution(tmp, weightName, runNumber, channelNumber);
- } else {
- Error("AddDistributionFromFile","%s is not a Histogram",histName.Data());
- throw std::runtime_error("Throwing 8");
- }
- }
- m_countingMode=false;
-
- delete rootFile;
-
- // Return to the original ROOT directory
- gDirectory = origDir;
-
- return 0;
-
-}
-Int_t Root::TPileupReweighting::AddDistributionFromFile(const TString fileName, const TString histName, Int_t runNumber, Int_t channelNumber) {
- return AddDistributionFromFile(fileName,histName,"pileup",runNumber,channelNumber);
-}
-
-Int_t Root::TPileupReweighting::AddLumiCalcFile(const TString fileName, const TString trigger) {
+Int_t CP::TPileupReweighting::AddLumiCalcFile(const TString fileName, const TString trigger) {
if(m_isInitialized) {
Error("AddLumiCalcFile","You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
@@ -831,44 +732,48 @@ Int_t Root::TPileupReweighting::AddLumiCalcFile(const TString fileName, const TS
m_lumicalcFiles[trigger] = fileName;
if(trigger=="None") {
Info("AddLumiCalcFile","Adding LumiMetaData (scale factor=%f)...",m_dataScaleFactorX);
+ //structure expected is as given by iLumiCalc:
+ // RunNbr, AvergeInteractionPerXing, IntLumi
+ UInt_t runNbr=0;Float_t intLumi=0;UInt_t lbn=0;TBranch *b_runNbr=0;TBranch *b_intLumi=0;TBranch *b_lbn=0;
+ Float_t mu=0.; TBranch *b_mu=0;
+ if(tmp->SetBranchAddress("RunNbr",&runNbr,&b_runNbr)!=0) {
+ Error("AddLumiCalcFile","Could not find RunNbr branch in Data TTree");throw std::runtime_error("Could not find RunNbr branch in Data TTree");
+ }
+ if(tmp->SetBranchAddress("AvergeInteractionPerXing",&mu,&b_mu)!=0) {
+ Error("AddLumiCalcFile","Could not find AvergeInteractionPerXing branch in Data TTree");throw std::runtime_error("Could not find AvergeInteractionPerXing branch in Data TTree");
+ }
+ if(tmp->SetBranchAddress("IntLumi",&intLumi,&b_intLumi)!=0) {
+ Error("AddLumiCalcFile","Could not find IntLumi branch in Data TTree");throw std::runtime_error("Could not find IntLumi branch in Data TTree");
+ }
+ if(tmp->SetBranchAddress("LBStart",&lbn,&b_lbn)!=0) {
+ Error("AddLumiCalcFile","Could not find LBStart branch in Data TTree");throw std::runtime_error("Could not find LBStart branch in Data TTree");
+ }
+ long nEntries = tmp->GetEntries();
+
+ for(long i=0;i<nEntries;i++) {
+ b_runNbr->GetEntry(i);b_intLumi->GetEntry(i);b_mu->GetEntry(i);
+ runNbr += m_lumicalcRunNumberOffset;
+ //save the lumi by run and lbn
+ b_lbn->GetEntry(i);
+
+ Run& r = m_runs[runNbr];
+ //rescale the mu value
+ mu *= m_dataScaleFactorX;
+ if(trigger=="None") {r.lumiByLbn[lbn].first += intLumi; r.lumiByLbn[lbn].second = mu;}
+ //fill into input data histograms
+ //check if we need to create an empty histogram
+
+ if(r.inputHists.find(trigger) == r.inputHists.end()) {
+ r.inputHists[trigger] = CloneEmptyHistogram(runNbr,-1);
+ }
+ r.inputHists[trigger]->Fill(mu,intLumi);
+ }
+ m_countingMode=false;
+ m_lumiVectorIsLoaded=true;
} else {
Info("AddLumiCalcFile","Adding LumiMetaData for DataWeight (trigger=%s) (scale factor=%f)...",trigger.Data(),m_dataScaleFactorX);
}
- //structure expected is as given by iLumiCalc:
- // RunNbr, AvergeInteractionPerXing, IntLumi
- UInt_t runNbr=0;Float_t intLumi=0;UInt_t lbn=0;TBranch *b_runNbr=0;TBranch *b_intLumi=0;TBranch *b_lbn=0;
- Float_t mu=0.; TBranch *b_mu=0;
- if(tmp->SetBranchAddress("RunNbr",&runNbr,&b_runNbr)!=0) {
- Error("AddLumiCalcFile","Could not find RunNbr branch in Data TTree");throw std::runtime_error("Could not find RunNbr branch in Data TTree");
- }
- if(tmp->SetBranchAddress("AvergeInteractionPerXing",&mu,&b_mu)!=0) {
- Error("AddLumiCalcFile","Could not find AvergeInteractionPerXing branch in Data TTree");throw std::runtime_error("Could not find AvergeInteractionPerXing branch in Data TTree");
- }
- if(tmp->SetBranchAddress("IntLumi",&intLumi,&b_intLumi)!=0) {
- Error("AddLumiCalcFile","Could not find IntLumi branch in Data TTree");throw std::runtime_error("Could not find IntLumi branch in Data TTree");
- }
- if(tmp->SetBranchAddress("LBStart",&lbn,&b_lbn)!=0) {
- Error("AddLumiCalcFile","Could not find LBStart branch in Data TTree");throw std::runtime_error("Could not find LBStart branch in Data TTree");
- }
- long nEntries = tmp->GetEntries();
- //std::map<Int_t, TH1*>& histograms = m_inputDataHistograms["pileup"];
- std::map<Int_t, TH1*>& histograms = (trigger=="None") ? m_inputHistograms["pileup"][-1] : m_triggerInputHistograms[trigger];
- for(long i=0;i<nEntries;i++) {
- b_runNbr->GetEntry(i);b_intLumi->GetEntry(i);b_mu->GetEntry(i);
- //save the lumi by run and lbn
- b_lbn->GetEntry(i);
- m_lumiByRunAndLbn[runNbr][lbn] = intLumi;
- //rescale the mu value
- mu *= m_dataScaleFactorX;
- //fill into input data histograms
- //check if we need to create an empty histogram
- if(histograms.find(runNbr) == histograms.end()) {
- histograms[runNbr] = CloneEmptyHistogram("pileup",runNbr,-1);
- }
- histograms[runNbr]->Fill(mu,intLumi);
- }
- m_countingMode=false;
- m_lumiVectorIsLoaded=true;
+
}
}
@@ -880,7 +785,7 @@ Int_t Root::TPileupReweighting::AddLumiCalcFile(const TString fileName, const TS
return 0;
}
-Int_t Root::TPileupReweighting::AddConfigFile(const TString fileName) {
+Int_t CP::TPileupReweighting::AddConfigFile(const TString fileName) {
//open the file and look for config TTrees
//recognised TTrees are: ChannelMetaData, MCPileupReweighting, MCCustomReweighting, DataCustomReweighting
@@ -899,19 +804,37 @@ Int_t Root::TPileupReweighting::AddConfigFile(const TString fileName) {
throw std::runtime_error(toThrow);
} else {
//try to get the the known TTrees
- TTree *tmp = (TTree*)rootFile->Get( "MCPileupReweighting" );
+ TTree *tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"MCPileupReweighting").c_str() );
+ if(tmp) {
+ Info("AddConfigFile","Adding MCPileupReweighting...");
+ AddDistributionTree(tmp,rootFile);
+ m_countingMode=false;
+ }
+ tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"MCCustomReweighting").c_str() );
+ if(tmp) {
+ Info("AddConfigFile","Adding MCCustomReweighting...");
+ AddDistributionTree(tmp,rootFile);
+ m_countingMode=false;
+ }
+ tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"DataCustomReweighting").c_str() );
+ if(tmp) {
+ Info("AddConfigFile","Adding DataCustomReweighting...");
+ AddDistributionTree(tmp,rootFile);
+ m_countingMode=false;
+ }
+ tmp=0; tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/MCPileupReweighting").c_str() ); //allow trees in the PileupReweighting dir instead
if(tmp) {
Info("AddConfigFile","Adding MCPileupReweighting...");
AddDistributionTree(tmp,rootFile);
m_countingMode=false;
}
- tmp = 0;tmp = (TTree*)rootFile->Get( "MCCustomReweighting" );
+ tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/MCCustomReweighting").c_str() );
if(tmp) {
Info("AddConfigFile","Adding MCCustomReweighting...");
AddDistributionTree(tmp,rootFile);
m_countingMode=false;
}
- tmp = 0;tmp = (TTree*)rootFile->Get( "DataCustomReweighting" );
+ tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/DataCustomReweighting").c_str() );
if(tmp) {
Info("AddConfigFile","Adding DataCustomReweighting...");
AddDistributionTree(tmp,rootFile);
@@ -927,1101 +850,759 @@ Int_t Root::TPileupReweighting::AddConfigFile(const TString fileName) {
return 0;
}
-void Root::TPileupReweighting::AddDataDistribution(const TString& dataRootFileName,
- const TString& dataRootHistName,Int_t runNumber) {
- Warning("AddDataDistribution","Method Depricated. Please use either AddConfigFile or AddDistributionFromFile");
- if(dataRootHistName=="LumiMetaData") AddLumiCalcFile(dataRootFileName);
- else AddDistributionFromFile(dataRootFileName,dataRootHistName,"pileup",runNumber,-1);
-}
-
-void Root::TPileupReweighting::AddMCDistribution(const TString& mcRootFileName,
- const TString& mcRootHistName, int runNumber, int channelNumber) {
- Warning("AddMCDistribution","Method Depricated. Please use either AddConfigFile or AddDistributionFromFile");
- if(mcRootHistName=="MCPileupReweighting") AddConfigFile(mcRootFileName);
- else AddDistributionFromFile(mcRootFileName,mcRootHistName,"pileup",runNumber,channelNumber);
-}
-
-void Root::TPileupReweighting::AddCustomDataDistribution(const TString& customRootFileName,
- const TString& customRootHistName, TString customName, int runNumber) {
- Warning("AddCustomDataDistribution","Method Depricated. Please use either AddConfigFile or AddDistributionFromFile");
- if(customRootHistName=="DataCustomReweighting") AddConfigFile(customRootFileName);
- else AddDistributionFromFile(customRootFileName,customRootHistName,customName,runNumber,-1);
-}
+Int_t CP::TPileupReweighting::Initialize() {
+ //Need to calculate a period weight for each period
+ //Need a reweighting histogram for each period
+ //for merged mc run numbers, we must generate or modify the existing mc distributions
-void Root::TPileupReweighting::AddCustomMCDistribution(const TString& customRootFileName,
- const TString& customRootHistName, TString customName, int runNumber, int channelNumber) {
- Warning("AddCustomMCDistribution","Method Depricated. Please use either AddConfigFile or AddDistributionFromFile");
- if(customRootHistName=="MCCustomReweighting") AddConfigFile(customRootFileName);
- else AddDistributionFromFile(customRootFileName,customRootHistName,customName,runNumber,channelNumber);
-}
-//channelNumber=-1 is data
-Int_t Root::TPileupReweighting::FactorizeDistribution(TH1* hist, const TString weightName,
- Int_t channelNumber, Int_t periodNumber,bool includeInMCRun,bool includeInGlobal) {
+ //print out the period assignments if in debug mode
+ if(m_debugging) {
+ PrintPeriods();
+ }
- //determine which mc runNumber this periodNumber is assigned to
- Int_t mcRunNumber = m_periodToMCRun[periodNumber];
- if(includeInGlobal) {
- globalTotals[weightName][channelNumber] += hist->GetSumOfWeights();
- globalNumberOfEntries[weightName][channelNumber] += hist->GetEntries();
- }
- if(includeInMCRun) {
- periodTotals[weightName][channelNumber][mcRunNumber] += hist->GetSumOfWeights();
- //initialize the corresponding data value if necessary
- if(periodTotals[weightName][-1].find(mcRunNumber)==periodTotals[weightName][-1].end()) {
- periodTotals[weightName][-1][mcRunNumber]=0;
+ if(m_countingMode) {
+ Info("Initialize","In Config File Generating mode. Remember to call WriteToFile!");
+ //need to check no periods have subperiods, this is not allowed in counting mode
+ for(auto period : m_periods) {
+ if(period.second->subPeriods.size()!=0) {
+ Error("Initialize","You cannot have subperiods when in Config File Generating mode");
+ throw std::runtime_error("Throwing 44: In Config File Generating mode, but detected subperiods in the period definition. This is not allowed.");
+ }
}
+ //histograms are instantiated in the event loop as the channels occur
+ //can delay period definition here
+ //but since we are here, check that the user has at least defined some periods .. i.e.
+// if(m_periods.size()==1) {
+// Error("Initialize", "In Config File Generating mode, but no periods have been defined. This makes no sense!? Define some periods please (with AddPeriod method)!");
+// throw std::runtime_error("Throwing 43: In Config File Generating mode, but no periods have been defined. This makes no sense!? Define some periods please (with AddPeriod method)!");
+// }
+ m_isInitialized=true;
+ return 0;
}
- periodTotals[weightName][channelNumber][periodNumber] += hist->GetSumOfWeights();
-
- //initialize the corresponding data value if necessary
- if(periodTotals[weightName][-1].find(periodNumber)==periodTotals[weightName][-1].end()) {
- periodTotals[weightName][-1][periodNumber]=0;
- }
-
- //add to the primary/secondary/tertiary histograms
- if(hist->InheritsFrom("TH3")) {
- //not supported yet
- Error("initialize","3d not supported yet");throw 40;
- } else if(hist->InheritsFrom("TH2")) {
- TH2* histo2d = dynamic_cast<TH2*>(hist);
- //project out and add that to primary
- TH1D* proj = histo2d->ProjectionX();
- if(!primaryDistributions[weightName][channelNumber][periodNumber]) {
- //need to make a copy
- TH1D* projCopy = dynamic_cast<TH1D*>(proj->Clone(proj->GetName()));
- projCopy->SetDirectory(0);
- primaryDistributions[weightName][channelNumber][periodNumber] = projCopy;
- //initialize the corresponding data histogram if necessary
- if(!primaryDistributions[weightName][-1][periodNumber]) {
- TH1D* projEmpty = dynamic_cast<TH1D*>(proj->Clone(proj->GetName()));
- projEmpty->Reset();
- projEmpty->SetDirectory(0);
- primaryDistributions[weightName][-1][periodNumber]=projEmpty;
+
+
+ //go through all periods, use inputHists to populate the primary and secondary dists, and fill sumOfWeights and numberOfEntries
+ for(auto period : m_periods) {
+ if(period.first != period.second->id) continue; //skips remappings
+ for(auto& inputHist : period.second->inputHists) {
+ int channel = inputHist.first;
+ TH1* hist = inputHist.second;
+ period.second->sumOfWeights[channel] += hist->GetSumOfWeights();
+ period.second->numberOfEntries[channel] += hist->GetEntries();
+ if(hist->GetDimension()==1) {
+ period.second->primaryHists[channel] = dynamic_cast<TH1D*>(CloneEmptyHistogram(period.first,channel));
+ period.second->primaryHists[channel]->Add(hist);
}
- } else {
- primaryDistributions[weightName][channelNumber][periodNumber]->Add(proj);
}
- if(!secondaryDistributions[weightName][channelNumber][periodNumber]) {
- secondaryDistributions[weightName][channelNumber][periodNumber] = dynamic_cast<TH2D*>(CloneEmptyHistogram(weightName,periodNumber,channelNumber));
- //initialize the corresponding data histogram if necessary
- if(!secondaryDistributions[weightName][-1][periodNumber]) {
- secondaryDistributions[weightName][-1][periodNumber]=dynamic_cast<TH2D*>(CloneEmptyHistogram(weightName,periodNumber,-1));
- }
- }
- secondaryDistributions[weightName][channelNumber][periodNumber]->Add(hist);
-
- if(includeInMCRun) {
- if(!primaryDistributions[weightName][channelNumber][mcRunNumber]) {
- //need to make a copy for the run number, to avoid double deleting
- TH1D* projCopy = dynamic_cast<TH1D*>(proj->Clone(proj->GetName()));
- projCopy->SetDirectory(0);
- primaryDistributions[weightName][channelNumber][mcRunNumber] = projCopy;
- //initialize the corresponding data histogram if necessary
- if(!primaryDistributions[weightName][-1][mcRunNumber]) {
- TH1D* projEmpty = dynamic_cast<TH1D*>(proj->Clone(proj->GetName()));
- projEmpty->Reset();
- projEmpty->SetDirectory(0);
- primaryDistributions[weightName][-1][mcRunNumber]=projEmpty;
- }
- } else {
- primaryDistributions[weightName][channelNumber][mcRunNumber]->Add(proj);
- }
- if(!secondaryDistributions[weightName][channelNumber][mcRunNumber]) {
- secondaryDistributions[weightName][channelNumber][mcRunNumber] = dynamic_cast<TH2D*>(CloneEmptyHistogram(weightName,mcRunNumber,channelNumber));
- //initialize the corresponding data histogram if necessary
- if(!secondaryDistributions[weightName][-1][mcRunNumber]) {
- secondaryDistributions[weightName][-1][mcRunNumber]=dynamic_cast<TH2D*>(CloneEmptyHistogram(weightName,mcRunNumber,-1));
- }
- }
- secondaryDistributions[weightName][channelNumber][mcRunNumber]->Add(hist);
- }
- delete proj;
-
- } else {
- //regular 1D distribution
- if(!primaryDistributions[weightName][channelNumber][periodNumber]) {
- primaryDistributions[weightName][channelNumber][periodNumber] = dynamic_cast<TH1D*>(CloneEmptyHistogram(weightName,periodNumber,channelNumber));
- //initialize the corresponding data histogram if necessary
- if(!primaryDistributions[weightName][-1][periodNumber]) {
- primaryDistributions[weightName][-1][periodNumber]=dynamic_cast<TH1D*>(CloneEmptyHistogram(weightName,periodNumber,-1));
- }
- }
- primaryDistributions[weightName][channelNumber][periodNumber]->Add(hist);
- if(includeInMCRun) {
- if(!primaryDistributions[weightName][channelNumber][mcRunNumber]) {
- primaryDistributions[weightName][channelNumber][mcRunNumber] = dynamic_cast<TH1D*>(CloneEmptyHistogram(weightName,mcRunNumber,channelNumber));
- //initialize the corresponding data histogram if necessary
- if(!primaryDistributions[weightName][-1][mcRunNumber]) {
- primaryDistributions[weightName][-1][mcRunNumber]=dynamic_cast<TH1D*>(CloneEmptyHistogram(weightName,mcRunNumber,-1));
- }
- }
- primaryDistributions[weightName][channelNumber][mcRunNumber]->Add(hist);
- }
- }
-
- return 0;
-}
-
-//=============================================================================
-// Initialization
-//=============================================================================
-Int_t Root::TPileupReweighting::initialize( const TString dataRootFileName,
- const TString dataRootHistName,
- const TString mcRootFileName,
- const TString mcRootHistName, int runNumber, int channel )
-{
- if(dataRootFileName != "") {
- AddMCDistribution(mcRootFileName,mcRootHistName,runNumber,channel);
- AddDataDistribution(dataRootFileName,dataRootHistName,runNumber);
- if(m_periods.size()==0) {
- Info("initialize","No Periods have been defined. I'm guessing you are doing old-style reweighting");
- Info("initialize","Creating dummy period and ignoring unrepresented data issues");
- AddPeriod(0,0,0); //assign mcRunNumber 0 to data runNumber 0-0.
- SetUnrepresentedDataAction(2);
- }
- }
- return Initialize();
-}
-
-Int_t Root::TPileupReweighting::Initialize() {
-
- //Need to calculate a period weight for each period
- //Need a reweighting histogram for each period
- //for merged mc run numbers, we must generate or modify the existing mc distributions
-
-
- //1. Deal with the registered mc mergers
- std::map<Int_t, bool> newMCRuns; //records if new mcRunNumbers have been created by merging
- for(std::map<Int_t, std::map<Int_t, Int_t> >::iterator mergers = m_mcMerges.begin();mergers != m_mcMerges.end();++mergers) {
- for(std::map<Int_t,Int_t>::iterator mergePair = (mergers->second).begin();mergePair!=(mergers->second).end();++mergePair) {
- Int_t run1 = mergePair->first;
- Int_t run2 = mergePair->second;
- Int_t runOut = mergers->first;
- if(run1==runOut) {
- if(m_debugging) Info("Initialize","Merging %d into %d",run2,run1);
- //merging run2 directly in to run1.
- //just combine run2 and run1 mc histos, and clear empty the run2 histograms
- for(std::map<TString, std::map<Int_t, std::map<Int_t,TH1*> > >::iterator mcHistos = m_inputHistograms.begin();mcHistos!=m_inputHistograms.end();++mcHistos) {
- TString weightName = mcHistos->first;
- for(std::map<Int_t, std::map<Int_t,TH1*> >::iterator channels = mcHistos->second.begin();channels!=mcHistos->second.end();++channels) {
- if((channels->first)<0) continue; //skip data channels
- if(channels->second.find(run2)!=channels->second.end()) {
- if(channels->second.find(run1)!=channels->second.end()) {
- channels->second[run1]->Add(channels->second[run2]);
- } else {
- channels->second[run1] = CloneEmptyHistogram(weightName,run1,channels->first);
- channels->second[run1]->Add(channels->second[run2]);
+ //also add all runs that belong to us ... only do the "None" trigger
+ for(auto& run : m_runs) {
+ if(!period.second->contains(run.first)) continue;
+ if(run.second.inputHists.find("None") == run.second.inputHists.end()) continue;
+ //add to list of runNumbers for this period
+ period.second->runNumbers.push_back(run.first);
+ TH1* hist = run.second.inputHists["None"];
+
+ //sweep over bins, checking for unrepresented data and acting accordingly
+ Int_t bin,binx,biny;
+ for(biny=1; biny<=hist->GetNbinsY(); biny++) {
+ for(binx=1; binx<=hist->GetNbinsX(); binx++) {
+ bin = hist->GetBin(binx,biny);
+ Double_t value = hist->GetBinContent(bin);
+ if(value==0.) continue;
+ //loop over channels, check for zero in this bin ... only need to do if we haven't checked ourselves
+ if(period.second->inputBinRedirect.find(bin) == period.second->inputBinRedirect.end()) {
+ for(auto& inputHist : period.second->inputHists) {
+ if(inputHist.first<0) continue; //skips data
+ Double_t mcValue = (inputHist.second)->GetBinContent(bin);
+ if(mcValue==0.) {
+ if(m_unrepresentedDataAction!=0 && m_debugging) Info("Initialize","Unrepresented data at coords [%f,%f] caused by periodNumber %d in channel %d",hist->GetXaxis()->GetBinCenter(binx),hist->GetYaxis()->GetBinCenter(biny),period.first,inputHist.first);
+ //if we are doing unrepaction=2, just set to 'not the bin number'. if we are doing action=3, find the nearest good bin
+ if(m_unrepresentedDataAction==0) {period.second->inputBinRedirect[bin] = bin-1;Error("Initialize","Unrepresented data at coords [%f,%f] caused by periodNumber %d in channel %d",hist->GetXaxis()->GetBinCenter(binx),hist->GetYaxis()->GetBinCenter(biny),period.first,inputHist.first);}
+ if(m_unrepresentedDataAction==1||m_unrepresentedDataAction==2) period.second->inputBinRedirect[bin] = bin-1;
+ else if(m_unrepresentedDataAction==3) period.second->inputBinRedirect[bin] = GetNearestGoodBin(period.first, bin);
+ break;
}
- channels->second[run2]->Reset();
}
+ if(period.second->inputBinRedirect.find(bin) == period.second->inputBinRedirect.end()) period.second->inputBinRedirect[bin] = bin; //a good bin
}
- }
- //merge the periods that were assigned
- //combine the periods
- if(m_periods.find(run2)!=m_periods.end()) {
- for(std::map<Int_t,std::pair<UInt_t,UInt_t> >::iterator periods2 = m_periods[run2].begin();periods2!=m_periods[run2].end();++periods2) {
- AddPeriod(run1,periods2->second.first,periods2->second.second);
- }
- m_periods.erase(run2);
- }
- //and remap run2 in to run1
- m_mcRemappings[run2]=run1;
- } else {
- //merging run1 and run2 in to a new run (run3).
- newMCRuns[runOut]=true;
- for(std::map<TString, std::map<Int_t, std::map<Int_t,TH1*> > >::iterator mcHistos = m_inputHistograms.begin();mcHistos!=m_inputHistograms.end();++mcHistos) {
- TString weightName = mcHistos->first;
- for(std::map<Int_t, std::map<Int_t,TH1*> >::iterator channels = mcHistos->second.begin();channels!=mcHistos->second.end();++channels) {
- if((channels->first)<0) continue;
- if(channels->second.find(run2)!=channels->second.end()) {
- if(channels->second.find(runOut)!=channels->second.end()) {
- channels->second[runOut]->Add(channels->second[run2]);
- } else {
- channels->second[runOut] = CloneEmptyHistogram(weightName,runOut,channels->first);
- channels->second[runOut]->Add(channels->second[run2]);
+ //if is a bad bin, act accordingly
+ if(period.second->inputBinRedirect[bin] != bin) {
+ run.second.badBins[bin]=value;
+ if(m_unrepresentedDataAction==1) {
+ //remove this data. we have to remember to do this for the trigger data too if this is the pileup weight
+ hist->SetBinContent(bin,0);
+ //also modify the m_inputTriggerHistograms ...
+ for(auto& triggerDists : run.second.inputHists) triggerDists.second->SetBinContent(bin,0);
+ } else if(m_unrepresentedDataAction==3) {
+ //reassign the content to that bin
+ hist->SetBinContent(period.second->inputBinRedirect[bin],hist->GetBinContent(period.second->inputBinRedirect[bin])+hist->GetBinContent(bin));
+ hist->SetBinContent(bin,0);
+ //also modify the m_inputTriggerHistograms ...
+ for(auto& triggerDists : run.second.inputHists) {
+ triggerDists.second->SetBinContent(period.second->inputBinRedirect[bin],triggerDists.second->GetBinContent(period.second->inputBinRedirect[bin])+triggerDists.second->GetBinContent(bin));
+ triggerDists.second->SetBinContent(bin,0);
}
- }
- if(channels->second.find(run1)!=channels->second.end()) {
- if(channels->second.find(runOut)!=channels->second.end()) {
- channels->second[runOut]->Add(channels->second[run1]);
- } else {
- channels->second[runOut] = CloneEmptyHistogram(weightName,runOut,channels->first);
- channels->second[runOut]->Add(channels->second[run1]);
- }
- }
- }
- }
- //merge the periods that were assigned, but do not clear the period assignment
- if(m_periods.find(run2)!=m_periods.end()) {
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods2 = m_periods[run2].begin();periods2!=m_periods[run2].end();++periods2) {
- AddPeriod(runOut,periods2->second.first,periods2->second.second);
- }
+ }
+ }
}
- if(m_periods.find(run1)!=m_periods.end()) {
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods2 = m_periods[run1].begin();periods2!=m_periods[run1].end();++periods2) {
- AddPeriod(runOut,periods2->second.first,periods2->second.second);
+ }
+ //create the period's 'data' hist if necessary
+ if( hist->GetDimension()==1 ) {
+ if(!period.second->primaryHists[-1] ) {
+ period.second->primaryHists[-1] = dynamic_cast<TH1D*>(CloneEmptyHistogram(period.first,-1));
+ period.second->triggerHists["None"] = dynamic_cast<TH1D*>(CloneEmptyHistogram(period.first,-1)); //this will remain unnormalized, unlike the 'primaryHist'
}
+ }
+ else if( hist->GetDimension()==2 ) {
+ if(!period.second->secondaryHists[-1] ) {
+ period.second->secondaryHists[-1] = dynamic_cast<TH2D*>(CloneEmptyHistogram(period.first,-1));
+ period.second->primaryHists[-1] = period.second->secondaryHists[-1]->ProjectionX(); period.second->primaryHists[-1]->SetDirectory(0);
+ period.second->primaryHists[-1]->Reset();
+ period.second->triggerHists["None"] = static_cast<TH1D*>(period.second->primaryHists[-1]->Clone("triggerHist"));
+ period.second->triggerHists["None"]->SetDirectory(0);
}
+
}
- }
- }
- //print out the period assignments if in debug mode
- if(m_debugging) {
- for(std::map<Int_t,std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator mcruns = m_periods.begin();mcruns!=m_periods.end();++mcruns) {
- Info("Initialize","MCRunNumber %d has following periods:",mcruns->first);
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods = mcruns->second.begin();periods!=mcruns->second.end();++periods) {
- Info("Initialize"," PeriodNumber %d: %d - %d",periods->first,periods->second.first,periods->second.second);
- }
}
+
}
- //Factorize the MC distributions to build up the totals and primary/seconday/tertiary distributions
- for(std::map<TString, std::map<Int_t, std::map<Int_t,TH1*> > >::iterator mcHistos = m_inputHistograms.begin();mcHistos!=m_inputHistograms.end();++mcHistos) {
- TString weightName = mcHistos->first;
- for(std::map<Int_t, std::map<Int_t,TH1*> >::iterator channels = mcHistos->second.begin();channels!=mcHistos->second.end();++channels) {
- Int_t channel = channels->first;
- if(channel<0) continue;
- for(std::map<Int_t,TH1*>::iterator mcruns = channels->second.begin();mcruns!=channels->second.end();++mcruns) {
- Int_t mcRunNum = mcruns->first;
- //all the periods associated to the mcrun will get the same distribution
- //only include the first subperiod in the totals and mcrunnum's distributions
- //and only include genuine mcrunnums (not ones generated in mergers) in global totals
- bool isFirst = true;
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator period = m_periods[mcRunNum].begin(); period!=m_periods[mcRunNum].end();++period) {
- FactorizeDistribution(mcruns->second,weightName,channel,period->first,isFirst,isFirst&&(!newMCRuns[mcRunNum]));
- isFirst=false;
+ Double_t unrepData(0);
+ //count up the unrepresented data, indexed by channel
+
+ for(auto& run : m_runs) {
+ for(auto& b : run.second.badBins) {
+ if(!b.second) continue;
+ int bin = b.first;
+ unrepData += b.second;
+ //find out which channel in which period caused this bad bin!
+ std::map<Int_t, bool> doneChannels;
+ for(auto& period : m_periods) {
+ if(period.first != period.second->id) continue;
+ if(!period.second->contains(run.first)) continue;
+ for(auto& inHist : period.second->inputHists) {
+ if(inHist.first<0) continue; //skips any data hists (shouldn't happen really)
+ if((inHist.second)->GetBinContent(bin)==0) { period.second->unrepData[inHist.first] += b.second; } //store unrep data by period too, because will use in GetUnrepresentedFraction
+ if(doneChannels[inHist.first]) continue; //dont doublecount the data if the channel was to blame across multiple periods
+ if((inHist.second)->GetBinContent(bin)==0) {unrepDataByChannel[inHist.first] += b.second;doneChannels[inHist.first]=true;}
}
}
}
}
- //build any trigger "data weight" histograms. Have one for each period, for each trigger
- for(std::map<Int_t, std::map<Int_t,std::pair<UInt_t,UInt_t> > >::iterator runPeriods = m_periods.begin();runPeriods != m_periods.end();++runPeriods) {
- for(std::map<Int_t,std::pair<UInt_t,UInt_t> >::iterator period = runPeriods->second.begin();period!=runPeriods->second.end();++period) {
- Int_t periodNumber = period->first;
- UInt_t pStart = period->second.first; UInt_t pEnd = period->second.second;
- if(m_inputHistograms.find("pileup")==m_inputHistograms.end()) continue; //need the "None" unprescaled mu distribution
- std::map<Int_t,TH1*> unprescaledDataHists = m_inputHistograms["pileup"][-1];
- m_triggerPrimaryDistributions["None"][periodNumber] = dynamic_cast<TH1D*>(CloneEmptyHistogram("pileup",periodNumber,-1));
- //loop over runs and see if we have corresponding in the unprescaled data map
- for(UInt_t r=pStart;r<=pEnd;r++) {
- if(unprescaledDataHists.find(r)!=unprescaledDataHists.end()) {
- if(unprescaledDataHists[r]) m_triggerPrimaryDistributions["None"][periodNumber]->Add(unprescaledDataHists[r]);
- }
- //also now look for trigger input histograms
- for(std::map<TString, std::map<Int_t, TH1*> >::iterator triggerDists = m_triggerInputHistograms.begin();triggerDists!=m_triggerInputHistograms.end();++triggerDists) {
- if(m_triggerPrimaryDistributions[triggerDists->first][periodNumber]==0) {
- m_triggerPrimaryDistributions[triggerDists->first][periodNumber] = dynamic_cast<TH1D*>(CloneEmptyHistogram("pileup",periodNumber,-1));
- Info("Initialize","Created Data Weight Histogram for [%s,%d]",triggerDists->first.Data(),periodNumber);
- }
- if(triggerDists->second.find(r)!=triggerDists->second.end()) {
- if(triggerDists->second[r]) {
- m_triggerPrimaryDistributions[triggerDists->first][periodNumber]->Add(triggerDists->second[r]);
- delete triggerDists->second[r]; triggerDists->second[r]=0;
- }
- }
- }
+ //now that all periods and runs have been checked, with redirects set up where necessary, we actually accumulate the data across the runs
+ //we should also check if there are any data runs that are not covered by the period assignments
+
+ double ignoredData(0);
+ for(auto& run : m_runs) {
+ bool covered(false);
+ TH1* hist = run.second.inputHists["None"];
+ for(auto& period : m_periods) {
+ if(!period.second->contains(run.first)) continue;
+ if(period.first!=-1) covered=true; //don't count the global period
+ if(hist->GetDimension()==1) {
+ period.second->primaryHists[-1]->Add(hist);
+ period.second->triggerHists["None"]->Add(hist);
+ } else if(hist->GetDimension()==2) {
+ period.second->secondaryHists[-1]->Add(hist);
+ TH1* proj = ((TH2*)hist)->ProjectionX();
+ period.second->primaryHists[-1]->Add(proj);
+ period.second->triggerHists["None"]->Add(proj);
+ delete proj;
}
+ period.second->sumOfWeights[-1] += hist->GetSumOfWeights();
+ period.second->numberOfEntries[-1] += hist->GetEntries();
}
+ if(!covered && m_periods.size()>1) {
+ Warning("Initialize","loaded data in run %d that is not covered by period assignments",run.first);
+ ignoredData += hist->GetSumOfWeights();
+ }
}
+ double totalData = (m_unrepresentedDataAction==1) ? (m_periods[-1]->sumOfWeights[-1]+unrepData) : m_periods[-1]->sumOfWeights[-1];
- std::map<TString,Double_t> unrepData;
- std::map<TString,std::map<Int_t,Bool_t> > doneRun; //keeps record of which data runs have been included in the global totals
-
- for(std::map<Int_t, std::map<Int_t,std::pair<UInt_t,UInt_t> > >::iterator runPeriods = m_periods.begin();runPeriods != m_periods.end();++runPeriods) {
- Int_t thisMCRunNumber = runPeriods->first;
- for(std::map<Int_t,std::pair<UInt_t,UInt_t> >::iterator period = runPeriods->second.begin();period!=runPeriods->second.end();++period) {
- Int_t periodNumber = period->first;
- UInt_t pStart = period->second.first; UInt_t pEnd = period->second.second;
- //for(std::map<TString,std::map<Int_t,TH1*> >::iterator dataHists = m_inputDataHistograms.begin();dataHists!=m_inputDataHistograms.end();++dataHists) {
- for(std::map<TString,std::map<Int_t,std::map<Int_t,TH1*> > >::iterator dataHists = m_inputHistograms.begin();dataHists!=m_inputHistograms.end();++dataHists) {
- TString weightName = dataHists->first;
- //get list of loaded runs associated to this period
- std::vector<UInt_t> associatedRuns;
- //std::map<Int_t,TH1*> hists = m_inputDataHistograms[weightName];
- std::map<Int_t,TH1*> hists = (dataHists->second)[-1];
- //look for loaded runs
- for(UInt_t run=pStart;run<=pEnd;run++) {
- if(hists.find(run)!=hists.end() && hists[run]) {
- associatedRuns.push_back(run);
- }
- }
- //loop over data runs in this period and build up the data distributions
- for(std::vector<UInt_t>::iterator run = associatedRuns.begin();run!=associatedRuns.end();++run) {
- TH1* thisHist = hists[*run];
- //look for unrepresented data in any of the bins of any of the corresponding channels
- for(std::map<Int_t,std::map<Int_t,TH1*> >::iterator channels=m_inputHistograms[weightName].begin();channels!=m_inputHistograms[weightName].end();++channels) {
- if((channels->first)<0) continue; //skip data channels
- if((channels->second)[thisMCRunNumber]) {
- Int_t bin,binx,biny,binz;
- for(binz=1; binz<=thisHist->GetNbinsZ(); binz++) {
- for(biny=1; biny<=thisHist->GetNbinsY(); biny++) {
- for(binx=1; binx<=thisHist->GetNbinsX(); binx++) {
- bin = thisHist->GetBin(binx,biny,binz);
- Double_t value = thisHist->GetBinContent(bin);
- if(value==0.) continue;
- Double_t mcValue = (channels->second)[thisMCRunNumber]->GetBinContent(bin);
- if(mcValue==0. && !m_badbins[weightName][*run][bin]) { //we don't want to double-count as we loop over channels
- //unrepresented data bin
- m_badbins[weightName][*run][bin]=true;
- unrepData[weightName] += value;
- if(m_unrepresentedDataAction==0 || m_debugging) {
- if(m_debugging) Error("Initialize","Unrepresented data at coords [%f,%f,%f] caused by mcrun %d in channel %d",thisHist->GetXaxis()->GetBinCenter(binx),thisHist->GetYaxis()->GetBinCenter(biny),thisHist->GetZaxis()->GetBinCenter(binz),thisMCRunNumber,channels->first);
- }
- if(m_unrepresentedDataAction==1) thisHist->SetBinContent(bin,0);
- }
- }
- }
- }
- }
- }
- dataPeriodRunTotals[weightName][periodNumber][*run] += thisHist->GetSumOfWeights(); //for random runnumber
- dataPeriodRunTotals[weightName][thisMCRunNumber][*run] += thisHist->GetSumOfWeights(); //for random runnumber
- //also record the binned totals - wont this use up loads of memory!!?? Do only for "pileup" weight then
- if(weightName=="pileup") {
- for(int binNum=0;binNum<=thisHist->GetNbinsX();++binNum) {
- if(thisHist->GetBinContent(binNum)==0) continue;
- dataPeriodBinnedRunTotals["pileup"][binNum][periodNumber][*run]+=thisHist->GetBinContent(binNum);
- dataPeriodBinnedRunTotals["pileup"][binNum][thisMCRunNumber][*run]+=thisHist->GetBinContent(binNum);
- //keep track of totals in the run=0 slot - used by GetRandomRunNumber
- dataPeriodBinnedRunTotals["pileup"][binNum][periodNumber][0]+=thisHist->GetBinContent(binNum);
- dataPeriodBinnedRunTotals["pileup"][binNum][thisMCRunNumber][0]+=thisHist->GetBinContent(binNum);
- }
- }
- //don't double-count the runs. we track the run totals in dataPeriodRunTotals["pileup"][-1]
- if(weightName=="pileup" && dataPeriodRunTotals["pileup"][-1].find(*run)==dataPeriodRunTotals["pileup"][-1].end()) {
- dataPeriodRunTotals["pileup"][-1][*run]=thisHist->GetSumOfWeights();
- int lumiVectorPosition = -1;
- int runNbr = *run;
- if(runNbr>=177986&&runNbr<=178109) lumiVectorPosition=0;//period B
- else if(runNbr>=179710&&runNbr<=180481) lumiVectorPosition=1;//period D
- else if(runNbr>=180614&&runNbr<=180776) lumiVectorPosition=2;//period E
- else if(runNbr>=182013&&runNbr<=182519) lumiVectorPosition=3;//period F
- else if(runNbr>=182726&&runNbr<=183462) lumiVectorPosition=4;//period G
- else if(runNbr>=183544&&runNbr<=184169) lumiVectorPosition=5;//period H
- else if(runNbr>=185353&&runNbr<=186493) lumiVectorPosition=6;//period I
- else if(runNbr>=186516&&runNbr<=186755) lumiVectorPosition=7;//period J
- else if(runNbr>=186873&&runNbr<=187815) lumiVectorPosition=8;//period K
- else if(runNbr>=188902&&runNbr<=190343) lumiVectorPosition=9; //period L
- else if(runNbr>=190503&&runNbr<=191933) lumiVectorPosition=10; //period M
- if(lumiVectorPosition>=0) { m_integratedLumiVector[lumiVectorPosition]+=(thisHist->GetSumOfWeights())/1E6; }
- }
- //we will manage the global total addition here... to prevent double counting when config file loaded overlapping periods
- //NOTE: May want to check at AddConfigFile stage if config file contains conflicting period config.... need to decide if it is ok or not!!??
- if(!doneRun[weightName][*run]) {
- globalTotals[weightName][-1] += thisHist->GetSumOfWeights();
- globalNumberOfEntries[weightName][-1] += thisHist->GetEntries();
- doneRun[weightName][*run] = true;
- }
- //factorize out the distribution and add it to existing distributions
- //include this distribution in the mcRunNumber distribution and in the global total (only if mcRunNumber is a genuine run num. Don't want to double-count)
- FactorizeDistribution(thisHist,weightName,-1,periodNumber,true,false/* !newMCRuns[thisMCRunNumber]*/);
- } //end loop over data in period
- } //loop over weights
- } //loop over periods in mcrun
- } //loop over mcruns
+ if(ignoredData>0.) Warning("Initialize", "Period Assignments missed %f%% data",100.*ignoredData/totalData);
- //we should also check if there are any data runs that are not covered by the period assignments
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > >::iterator weights=m_inputHistograms.begin();weights!=m_inputHistograms.end();++weights) {
- Double_t ignoredData = 0;
- if(weights->second.find(-1)==weights->second.end()) continue; //no data
- for(std::map<Int_t,TH1*>::iterator runs = weights->second[-1].begin();runs!=weights->second[-1].end();++runs) {
- //look for period assignment for run
- Int_t periodNum=-1;
- for(std::map<Int_t, std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator mcruns = m_periods.begin();mcruns!=m_periods.end();++mcruns) {
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=mcruns->second.begin();periods!=mcruns->second.end();++periods) {
- if((UInt_t)runs->first >= periods->second.first && (UInt_t)runs->first <= periods->second.second) {
- periodNum = periods->first;break;
- }
- }
- if(periodNum!=-1) break;
- }
- if(periodNum==-1) {
- //this run isn't assigned to a period. Give warning
- Warning("Initialize","%s weight contains loaded data in run %d that is not covered by period assignments",weights->first.Data(),runs->first);
- ignoredData += (runs->second)->GetSumOfWeights();
+
+
+ if(unrepData) {
+ double frac = unrepData/totalData;
+ if( frac > m_unrepDataTolerance) {
+ Error("Initialize", "%f%% unrepresented data, which suggests something is wrong with your prw config. Try EnableDebugging(true) to investigate",100.* (unrepData/m_periods[-1]->sumOfWeights[-1]));
+ }
+ if(m_unrepresentedDataAction==1) {
+ Warning("Initialize","has %f%% unrepresented data. This was removed (UnrepresentedDataAction=1)",100.*frac);
+ } else if(m_unrepresentedDataAction==2) {
+ Warning("Initialize","has %f%% unrepresented data. This was kept in (UnrepresentedDataAction=2)",100.*frac);
+ } else if(m_unrepresentedDataAction==3) {
+ Warning("Initialize","has %f%% unrepresented data. This was reassigned (UnrepresentedDataAction=3)",100.*frac);
+ } else if(m_unrepresentedDataAction==0) {
+ Error("Initialize","has %f%% unrepresented data:",100.*frac);
+ //print the report of which channels caused it
+ for(auto& it : unrepDataByChannel) {
+ Error("Initialize"," Channel %d caused %f%% of the unrepresented data (nEntries=%d)",it.first,100.*it.second/unrepData,m_periods[-1]->numberOfEntries[it.first]);
}
+ Error("Initialize","Exiting. You must decide how to proceed...");
+ Error("Initialize","1) use AddPeriod or RemapPeriod to redefine the mc periods to include this data. You should not need to regenerate your mc config file");
+ Error("Initialize","2) use SetUnrepresentedDataAction(1) to remove this data from the weight calculations. You should also veto such data events (using IsUnrepresentedData(..,..) method)");
+ Error("Initialize","3) use SetUnrepresentedDataAction(2) to leave this data in the calculation. I hope you know what you're doing!!");
+ Error("Initialize","4) use SetUnrepresentedDataAction(3) to reassign the data to the nearest representable bin");
+ throw std::runtime_error("Throwing exception 22: Unrepresented data exists. Make a choice how to handle this. See err log for details");
}
- if(ignoredData>0.) Warning("Initialize", "Period Assignments missed %f data in %s weight",ignoredData,weights->first.Data());
- }
-
- bool displayOptions = false;
- bool toleranceViolated = false;
- for(std::map<TString,Double_t>::iterator it=unrepData.begin();it!=unrepData.end();++it) {
- if(it->second>0.) {
- //check if unrepresented data violates tolerance criteria, regardless of action mode
- if((it->second/globalTotals[it->first][-1]) > m_unrepDataTolerance) {
- Error("Initialize", "%s weight has %f%% unrepresented data, which suggests something is wrong with your prw config. Try EnableDebugging(true) to investigate",it->first.Data(),100.*it->second/globalTotals[it->first][-1]);
- toleranceViolated=true;
- }
- if(m_unrepresentedDataAction==0 || m_debugging) {
- Error("Initialize","%s weight has %f%% unrepresented data",it->first.Data(),100.*it->second/globalTotals[it->first][-1]);
- //sweep back over the bad bins and add up the unrepresented data grouped by channel
- std::map<Int_t,Double_t> unrepChannelTotals;
- for(std::map<Int_t, std::map<Int_t, Bool_t> >::iterator badRuns=m_badbins[it->first].begin();badRuns!=m_badbins[it->first].end();++badRuns) {
- //find out which mcrunNum this run belongs to
- Int_t mcrunNum=-1;
- for(std::map<Int_t, std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator mcruns = m_periods.begin();mcruns!=m_periods.end();++mcruns) {
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::iterator periods=mcruns->second.begin();periods!=mcruns->second.end();++periods) {
- if((UInt_t)badRuns->first >= periods->second.first && (UInt_t)badRuns->first <= periods->second.second) {
- mcrunNum = mcruns->first;break;
- }
- }
- if(mcrunNum!=-1) break;
- }
- //now sweep over mc channels and check this periodNum for the bad bins
- for(std::map<Int_t, std::map<Int_t, TH1*> >::iterator channels=m_inputHistograms[it->first].begin();channels!=m_inputHistograms[it->first].end();++channels) {
- if(channels->first < 0) continue;
- for(std::map<Int_t, Bool_t>::iterator bins=badRuns->second.begin();bins!=badRuns->second.end();++bins) {
- if(bins->second && (channels->second)[mcrunNum]->GetBinContent(bins->first)==0.) {
- //found one of the bad bins.
- unrepChannelTotals[channels->first] += m_inputHistograms[it->first][-1][badRuns->first]->GetBinContent(bins->first);
- }
- }
- }
- }
- //display the results of the scan
- for(std::map<Int_t,Double_t>::iterator badChans = unrepChannelTotals.begin();badChans!=unrepChannelTotals.end();++badChans) {
- if(m_unrepresentedDataAction!=0 && (badChans->second*100./it->second)<0.000001) continue; //don't display the tiny channel errors if just debugging
- //quickly count up the number of entries in this channel
- long nEntries(0);
- for(std::map<Int_t,TH1*>::iterator cIt = m_inputHistograms[it->first][badChans->first].begin();cIt!=m_inputHistograms[it->first][badChans->first].end();++cIt) {
- nEntries += cIt->second->GetEntries();
- }
- Error("Initialize"," Channel %d has %f%% of the unrepresented (nEntries=%ld)",badChans->first,badChans->second*100./it->second,nEntries);
- }
- if(m_unrepresentedDataAction==0) displayOptions=true;
- } else if(m_unrepresentedDataAction==1) {
- Warning("Initialize","%s weight has %f%% unrepresented data. This was removed (UnrepresentedDataAction=1)",it->first.Data(),100.*it->second/(globalTotals[it->first][-1]+it->second));
- } else if(m_unrepresentedDataAction==2) {
- Warning("Initialize","%s weight has %f%% unrepresented data. This was kept in (UnrepresentedDataAction=2)",it->first.Data(),100.*it->second/(globalTotals[it->first][-1]));
- }
+ //if get here, through the exception if frac is bad
+ if( frac > m_unrepDataTolerance) {
+ throw std::runtime_error("Throwing exception 222: Some channel had too much unrepresented data. You should fix your prw file");
}
}
- if(toleranceViolated) {
- Error("Initialize","Some channel had over %f%% of unrepresented data. You should fix your prw file",m_unrepDataTolerance*100.);
- throw std::runtime_error("Throwing exception 222: Some channel had too much unrepresented data. You should fix your prw file");
- }
- if(displayOptions) {
- Error("Initialize","Exiting. You must decide how to proceed...");
- Error("Initialize","1) use AddPeriod or MergeMCRunNumbers to redefine the mc periods to include this data. You should not need to regenerate your mc config file");
- Error("Initialize","2) use SetUnrepresentedDataAction(1) to remove this data from the weight calculations. You should also veto such data events (using IsUnrepresentedData(..,..) method)");
- Error("Initialize","3) use SetUnrepresentedDataAction(2) to leave this data in the calculation. I hope you know what you're doing!!");
- throw std::runtime_error("Throwing exception 22: Unrepresented data exists. Make a choice how to handle this. See err log for details");
- }
-
- if(m_debugging) Info("Initialize","Normalizing 1D histograms...");
+ if(m_debugging) Info("Initialize","Normalizing histograms and cleaning up...");
//now that all the distributions are built. Normalize them all
- for(std::map<TString, std::map<Int_t, std::map<Int_t, TH1D*> > >::iterator weights=primaryDistributions.begin();weights!=primaryDistributions.end();++weights) {
- for(std::map<Int_t, std::map<Int_t, TH1D*> >::iterator channels=weights->second.begin();channels!=weights->second.end();++channels) {
- for(std::map<Int_t,TH1D*>::iterator mcrunsAndPeriods=channels->second.begin();mcrunsAndPeriods!=channels->second.end();++mcrunsAndPeriods) {
- normalizeHistogram(mcrunsAndPeriods->second);
- }
+ for(auto period : m_periods) {
+ if(period.first != period.second->id) continue;
+ for(auto pHist : period.second->primaryHists) {
+ normalizeHistogram(pHist.second);
}
- }
- if(m_debugging) Info("Initialize","Normalizing 2D histograms...");
- for(std::map<TString, std::map<Int_t, std::map<Int_t, TH2D*> > >::iterator weights=secondaryDistributions.begin();weights!=secondaryDistributions.end();++weights) {
- for(std::map<Int_t, std::map<Int_t, TH2D*> >::iterator channels=weights->second.begin();channels!=weights->second.end();++channels) {
- for(std::map<Int_t,TH2D*>::iterator mcrunsAndPeriods=channels->second.begin();mcrunsAndPeriods!=channels->second.end();++mcrunsAndPeriods) {
- normalizeHistogram(mcrunsAndPeriods->second);
- }
+ for(auto pHist : period.second->secondaryHists) {
+ normalizeHistogram(pHist.second);
}
- }
-
- if(m_debugging) Info("Initialize","Deleting input histograms...");
- //to save space, delete all the input histograms
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > >::iterator it0=m_inputHistograms.begin();it0!=m_inputHistograms.end();++it0) {
- for(std::map<Int_t,std::map<Int_t, TH1*> >::iterator it1=it0->second.begin();it1!=it0->second.end();++it1) {
- for(std::map<Int_t, TH1*>::iterator it2=it1->second.begin();it2!=it1->second.end();++it2) {
- delete it2->second;
- }
- it1->second.clear();
+ for(auto pHist : period.second->inputHists) {
+ delete pHist.second;
}
- it0->second.clear();
+ period.second->inputHists.clear();
}
- m_inputHistograms.clear();
+
+ //we keep the inputHists for data, because can use that to do random run numbers based on mu
+
if(m_debugging) Info("Initialize","...Done");
//if no input histograms were added, we are in counting mode
- if(m_countingMode) {
- Info("Initialize","In Config File Generating mode. Remember to call WriteToFile!");
- //histograms are instantiated in the event loop as the channels occur
- //but since we are here, check that the user has at least defined some periods
- if(m_periods.size()==0) {
- Error("Initialize", "In Config File Generating mode, but no periods have been defined. This makes no sense!? Define some periods please (with AddPeriod method)!");
- throw std::runtime_error("Throwing 43: In Config File Generating mode, but no periods have been defined. This makes no sense!? Define some periods please (with AddPeriod method)!");
- }
- }
+
m_isInitialized=true;
return 0;
}
-Bool_t Root::TPileupReweighting::IsUnrepresentedData(const TString weightName, Int_t runNumber, Float_t x, Float_t y, Float_t z) {
- if(m_emptyHistograms.find(weightName)==m_emptyHistograms.end()) {
- Error("IsUnrepresentedData", "Unknown weight %s", weightName.Data());
- throw std::runtime_error("Throwing 23: Unknown weight");
- }
- TH1* hist = m_emptyHistograms[weightName];
- Int_t bin = hist->GetXaxis()->FindFixBin(x);
- if(hist->GetDimension()==2) {
- Int_t biny = hist->GetYaxis()->FindFixBin(y);
- Int_t nx = hist->GetXaxis()->GetNbins()+2;
- bin += nx*biny;
- }
- if(hist->GetDimension()==3) {
- Int_t biny = hist->GetYaxis()->FindFixBin(y);
- Int_t nx = hist->GetXaxis()->GetNbins()+2;
- Int_t ny = hist->GetYaxis()->GetNbins()+2;
- Int_t binz = hist->GetZaxis()->FindFixBin(z);
- bin += nx*(biny+ny*binz);
- }
- return m_badbins[weightName][runNumber][bin];
-}
+//loop over the mc histograms, find bin closest to given bin that is "ok" (has all nonzero)
+//0 = ok bin, 1 = bad bin, 2 = out of range bin
+Int_t CP::TPileupReweighting::IsBadBin(Int_t periodNumber, Int_t bin) {
+
+ if(bin<0) return 2; //out of range
+
+ Period* p = m_periods[periodNumber];
+ for(auto& inHist : p->inputHists) {
+ int channelNumber = inHist.first;
+ if(channelNumber<0) continue; //skip data hists
+ if(bin > (inHist.second->GetNbinsX()+2)*(inHist.second->GetNbinsY()+2)*(inHist.second->GetNbinsZ()+2)) return 2; //definitely out of range by this point!
+ if(inHist.second->GetBinContent(bin)==0) return 1;
+ }
+ return 0;
-Bool_t Root::TPileupReweighting::IsUnrepresentedData(Int_t runNumber, Float_t x, Float_t y, Float_t z) {
- return IsUnrepresentedData("pileup",runNumber,x,y,z);
}
-Int_t Root::TPileupReweighting::WriteToFile(const TString fname) {
-
- if(!m_countingMode) {Warning("WriteToFile","Not in counting mode, so no file will be written");return 0;}
+//loop over the mc histograms, find bin closest to given bin that is "ok" (has all nonzero)
+Int_t CP::TPileupReweighting::GetNearestGoodBin(Int_t thisMCRunNumber, Int_t bin) {
+ int binDistance = 1;
+ do {
+ int res1 = IsBadBin(thisMCRunNumber,bin+binDistance);
+ if(!res1) return bin+binDistance;
+ int res2 = IsBadBin(thisMCRunNumber,bin-binDistance);
+ if(!res2) return bin-binDistance;
+ if(res1==2 && res2==2) { //run out of bins
+ Error("GetNearestGoodBin", "None of the bins are good!!??");
+ return -1;
+ }
+ binDistance++;
+ } while(true); //scary!
- //build a TTree with the correct mc structure, and dump the mc histogram info in to it
- //also build aggregate histograms across all channels - shouldn't be used as input histograms unless for a single channel
- TString filename = fname;
- filename += (filename=="") ? TString(this->GetName()) + ".prw.root" : "";
-
- //loop over the weights. "pileup" gets it own MCPileupReweighting ttree.
- //data goes in to DataCustomReweighting. other goes in to MCCustomReweighting
+ return -1;
- TFile* outFile = TFile::Open(filename,"RECREATE");
- Int_t r = WriteToFile(outFile);
- outFile->Close();
- return r;
}
-Int_t Root::TPileupReweighting::WriteToFile(TFile* outFile) {
- TDirectory* origDir = gDirectory;
-
- outFile->cd();
+Double_t CP::TPileupReweighting::GetUnrepresentedDataFraction(Int_t periodNumber,Int_t channel) {
- TTree *outTreeMC=0;
- TTree *outTreeCustomMC=0;
- TTree *outTreeData=0;
- TTree *outTreeCustomData=0;
- TTree *outTree = 0;//points to tree currently being written to
- Int_t channel = 0;UInt_t runNumber = 0;
- std::vector<UInt_t>* pStarts = 0;std::vector<UInt_t>* pEnds = 0;
- Char_t histName[150];Char_t customName[150];
+ //return unrepDataByChannel[channel]/GetSumOfEventWeights(-1); //the denominator is equivalent to asking for the total data!
- bool lastChannelWasData = true;
- for(std::map<TString, std::map<Int_t, std::map<Int_t, TH1*> > >::iterator weights=m_inputHistograms.begin();weights!=m_inputHistograms.end();++weights) {
- TString weightName = weights->first;
- strcpy(customName,weightName);
- if(outTree) {
- outTree=0;
- }
- for(std::map<Int_t, std::map<Int_t, TH1*> >::iterator channels=weights->second.begin();channels!=weights->second.end();++channels) {
- channel = channels->first;
- if(lastChannelWasData&&channel>=0) {
- if(outTree) {outTree=0;}
- lastChannelWasData=false;
- }
- if(channel<0) lastChannelWasData=true;
- if(weightName=="pileup") {
- if(channel>=0 && !outTree) {
- if(!outTreeMC) {
- outTreeMC = new TTree("MCPileupReweighting","MCPileupReweighting");
- outTreeMC->Branch("Channel",&channel);
- outTreeMC->Branch("RunNumber",&runNumber);
- outTreeMC->Branch("PeriodStarts",&pStarts);
- outTreeMC->Branch("PeriodEnds",&pEnds);
- outTreeMC->Branch("HistName",&histName,"HistName/C");
- }
- outTree = outTreeMC;
- } else if(channel<0 && !outTree) {
- if(!outTreeData) {
- outTreeData = new TTree("DataPileupReweighting","DataPileupReweighting");
- outTreeData->Branch("RunNumber",&runNumber);
- outTreeData->Branch("HistName",&histName,"HistName/C");
- }
- outTree = outTreeData;
- }
- } else {
- if(channel>=0 && !outTree) {
- if(!outTreeCustomMC) {
- outTreeCustomMC = new TTree("MCCustomReweighting","MCCustomReweighting");
- outTreeCustomMC->Branch("CustomName",&customName,"CustomName/C");
- outTreeCustomMC->Branch("Channel",&channel);
- outTreeCustomMC->Branch("RunNumber",&runNumber);
- outTreeCustomMC->Branch("PeriodStarts",&pStarts);
- outTreeCustomMC->Branch("PeriodEnds",&pEnds);
- outTreeCustomMC->Branch("HistName",&histName,"HistName/C");
- }
- outTree = outTreeCustomMC;
- } else if(channel<0 && !outTree) {
- if(!outTreeCustomData) {
- outTreeCustomData = new TTree("DataCustomReweighting","DataCustomReweighting");
- outTreeCustomData->Branch("CustomName",&customName,"CustomName/C");
- outTreeCustomData->Branch("RunNumber",&runNumber);
- outTreeCustomData->Branch("HistName",&histName,"HistName/C");
- }
- outTree = outTreeCustomData;
- }
- }
- for(std::map<Int_t, TH1*>::iterator runs=channels->second.begin();runs!=channels->second.end();++runs) {
- runNumber = runs->first;
- if(!lastChannelWasData) { //load periods only for mc
- pStarts = new std::vector<UInt_t>;pEnds = new std::vector<UInt_t>;
- for(std::map<Int_t,std::pair<UInt_t,UInt_t> >::iterator periods=m_periods[runNumber].begin();periods!=m_periods[runNumber].end();++periods) {
- pStarts->push_back(periods->second.first);pEnds->push_back(periods->second.second);
- }
- }
- strcpy(histName,runs->second->GetName());
- outTree->Fill();
- runs->second->Write();
- if(!lastChannelWasData) { //delete periods only for mc
- delete pStarts;delete pEnds;
- }
- }
- }
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetUnrepresentedDataFraction","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- //write the non-zero ttrees
- if(outTreeMC) {outTreeMC->Write();delete outTreeMC;outTreeMC=0;}
- if(outTreeData) {outTreeData->Write();delete outTreeData;outTreeData=0;}
- if(outTreeCustomMC) {outTreeCustomMC->Write();delete outTreeCustomMC;outTreeCustomMC=0;}
- if(outTreeCustomData) {outTreeCustomData->Write();delete outTreeCustomData;outTreeCustomData=0;}
+ //look for channel through sumOfWeights - because if the channel has no unrepData, it will have no entry in the unrepData map
+ if(p->sumOfWeights.find(channel) == p->sumOfWeights.end()) channel = m_periods[periodNumber]->defaultChannel;
- Info("WriteToFile", "Successfully generated config file: %s",outFile->GetName());
- Info("WriteToFile", "Happy Reweighting :-)");
- gDirectory = origDir;
+ if(p->sumOfWeights.find(channel) == p->sumOfWeights.end()) {
+ Error("GetUnrepresentedDataFraction","Unrecognised channel: %d",channel);
+ throw std::runtime_error("GetUnrepresentedDataFraction: Unrecognised channel");
+ }
- return 0;
-}
+ return p->unrepData[channel]/p->sumOfWeights[-1];
-Double_t Root::TPileupReweighting::GetDataWeight(Int_t runNumber, TString trigger) {
- //special mu-independent version of GetDataWeight. Will just use the global luminosity
- m_doGlobalDataWeight=true;
- double out = GetDataWeight(runNumber, trigger, 0);
- m_doGlobalDataWeight=false;
- return out;
}
-Double_t Root::TPileupReweighting::GetDataWeight(Int_t runNumber, TString trigger, Double_t x) {
- if(!m_isInitialized) { Error("GetDataWeight","Tool is not initialized"); return 0;}
- if(m_countingMode) return 0.;
- if(m_countingMode) {
- Error("GetDataWeight","You cannot get a weight when in config file generating mode. Please use FillWeight method");
- throw std::runtime_error("GetDataWeight: You cannot get a weight when in config file generating mode. Please use FillWeight method");
- }
- //determine which mc run number this is
- Int_t periodNumber = GetFirstFoundPeriodNumber(runNumber);
- if(periodNumber==-1) {
- Error("GetDataWeight","Could not find corresponding period for runNumber %d",runNumber);
- throw std::runtime_error("GetDataWeight: Could not find period for runNumber");
- }
- if(m_triggerPrimaryDistributions.find(trigger)==m_triggerPrimaryDistributions.end()) {
- //try to do a subtrigger calculation we need to reopen all the lumicalc files and do the fiddily prescaled luminosity calculation
- //will throw errors if can't
- CalculatePrescaledLuminosityHistograms(trigger);
+
+
+UInt_t CP::TPileupReweighting::GetRandomRunNumber(Int_t periodNumber) {
+ if(!m_isInitialized) { Info("GetRandomRunNumber","Initializing the subtool.."); Initialize(); }
+ if(m_countingMode) { return 0; } //do nothing when in counting mode
+
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetRandomRunNumber","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- TH1D* denomHist = m_triggerPrimaryDistributions[trigger][periodNumber];
- if(denomHist==0) {
- Error("GetDataWeight","Could not find trigger %s",trigger.Data());
- throw std::runtime_error("GetDataWeight: Could not find trigger");
+
+ double lumi = p->sumOfWeights[-1] * m_random3->Rndm();
+
+ //go through the runs assigned to this period ...
+ double lumiSum=0;
+ for(auto runNum : p->runNumbers) {
+ lumiSum += m_runs[runNum].inputHists["None"]->GetSumOfWeights();
+ if(lumiSum >= lumi) return runNum;
}
- TH1D* numerHist = m_triggerPrimaryDistributions["None"][periodNumber];
- if(numerHist==0) {
- Error("GetDataWeight","Could not find unprescaled trigger");
- throw std::runtime_error("GetDataWeight: Could not find unprescaled trigger. Please AddLumiCalc with a 'None' trigger");
+ Error("GetRandomRunNumber","overran integrated luminosity for periodNumber=%d",periodNumber);
+ throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomRunNumber");
+ return 0;
+}
+
+UInt_t CP::TPileupReweighting::GetRandomRunNumber(Int_t periodNumber, Double_t x) {
+ if(!m_isInitialized) { Info("GetRandomRunNumber","Initializing the subtool.."); Initialize(); }
+ if(m_countingMode) { return 0; } //do nothing when in counting mode
+
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetRandomRunNumber","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- if(m_doGlobalDataWeight) return numerHist->Integral(0,numerHist->GetNbinsX()+1)/denomHist->Integral(0,denomHist->GetNbinsX()+1);
+ int bin = m_emptyHistogram->FindFixBin(x);
- Int_t bin=numerHist->FindFixBin(x);
- return numerHist->GetBinContent(bin)/denomHist->GetBinContent(bin);
+ double lumi = p->primaryHists[-1]->GetBinContent(bin) * p->sumOfWeights[-1] * m_random3->Rndm();
-}
+ if(!lumi) { return 0; } //no lumi with that mu
-/** alternative random runnumber method for people who want to use custom reweighting config files to determine the random numbers */
-//Could tidy up code by removing the version below, and defaulting weightName to "pileup"
-UInt_t Root::TPileupReweighting::GetRandomRunNumber(Int_t mcRunNumber, TString weightName) {
- if(!m_isInitialized) { Error("GetRandomRunNumber","Tool is not initialized"); return 0;}
- if(m_countingMode) { return 0; } //do nothing when in counting mode
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
- const std::map<Int_t, Double_t>::const_iterator it = periodTotals[weightName][-1].find(mcRunNumber);
- if(it==periodTotals[weightName][-1].end()|| it->second == 0) {
- if(m_SetWarnings) Warning("GetRandomRunNumber","%d has no associated custom data. Returning 0",mcRunNumber);
- return 0; /*throw 33;*/
- }
- double lumi = it->second * m_random3->Rndm();
- double lumisum = 0.;
- std::map<UInt_t, Double_t>::const_iterator itend = dataPeriodRunTotals[weightName][mcRunNumber].end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = dataPeriodRunTotals[weightName][mcRunNumber].begin(); it2 != itend; ++it2)
- {
- lumisum += it2->second;
- if (lumisum >= lumi)
- return it2->first;
- }
- Error("GetRandomRunNumber","overran integrated luminosity for mcRunNumber=%d in custom weight=%s",mcRunNumber,weightName.Data());
- throw std::runtime_error("Throwing 46.1: overran integrated luminosity for mcRunNumber in custom weight");
- return 0;
-}
-/** generate random run number from lumicalc mu distributions */
-UInt_t Root::TPileupReweighting::GetRandomRunNumber(Int_t mcRunNumber) {
- if(!m_isInitialized) { Error("GetRandomRunNumber","Tool is not initialized"); return 0;}
- if(m_countingMode) { return 0; } //do nothing when in counting mode
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
- const std::map<Int_t, Double_t>::const_iterator it = periodTotals["pileup"][-1].find(mcRunNumber);
- if(it==periodTotals["pileup"][-1].end()|| it->second == 0) {
- if(m_SetWarnings) Warning("GetRandomRunNumber","%d has no associated lumi. Returning 0",mcRunNumber);
- return 0; /*throw 33;*/
- }
- double lumi = it->second * m_random3->Rndm();
- double lumisum = 0.;
- std::map<UInt_t, Double_t>::const_iterator itend = dataPeriodRunTotals["pileup"][mcRunNumber].end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = dataPeriodRunTotals["pileup"][mcRunNumber].begin(); it2 != itend; ++it2)
- {
- lumisum += it2->second;
- if (lumisum >= lumi)
- return it2->first;
- }
- Error("GetRandomRunNumber","overran integrated luminosity for mcRunNumber=%d",mcRunNumber);
- throw std::runtime_error("Throwing 46: overran integrated luminosity for mcRunNumber");
- return 0;
+ //go through the runs assigned to this period ...
+ double lumiSum=0;
+ for(auto runNum : p->runNumbers) {
+ lumiSum += m_runs[runNum].inputHists["None"]->GetBinContent(bin);
+ if(lumiSum >= lumi) return runNum;
+ }
+ Error("GetRandomRunNumber","overran integrated luminosity for periodNumber=%d",periodNumber);
+ throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomRunNumber");
+ return 0;
}
-///return random run numbers but only for runs where the mu is present. If the mu has no associated data, returns 0
-UInt_t Root::TPileupReweighting::GetRandomRunNumber(Int_t mcRunNumber, Double_t mu) {
- if(!m_isInitialized) { Error("GetRandomRunNumber","Tool is not initialized"); return 0;}
+UInt_t CP::TPileupReweighting::GetRandomLumiBlockNumber(UInt_t runNumber) {
+ if(!m_isInitialized) {Info("GetRandomLumiBlockNumber","Initializing the subtool.."); Initialize(); }
if(m_countingMode) { return 0; } //do nothing when in counting mode
- //determine the mubin
- TH1* hist = m_emptyHistograms["pileup"];
- if(!hist) { Error("GetRandomRunNumber", "Cannot do this without a lumicalc file!"); return 0; }
- Int_t muBin = hist->FindFixBin(mu);
-
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
+ double lumi = GetIntegratedLumi(runNumber,runNumber) * m_random3->Rndm() * 1E6 /* dont forget the lumi was divided by a million to get to pb */;
+ double lumisum = 0.;
- //check that this mcRunNumber actually has some lumi with this mu
- std::map<UInt_t, std::map<Int_t, std::map<UInt_t, Double_t> > >::iterator myMap1 = dataPeriodBinnedRunTotals["pileup"].find(muBin);
- if(myMap1 == dataPeriodBinnedRunTotals["pileup"].end()) {
- //no lumi with this mu
- return 0;
- }
- std::map<Int_t, std::map<UInt_t, Double_t> >::iterator myMap2 = myMap1->second.find(mcRunNumber);
- if(myMap2 == myMap1->second.end()) {
- //no lumi with this mu in this mcRunNumber
- return 0;
+ for(auto& lbn : m_runs[runNumber].lumiByLbn) {
+ if(m_unrepresentedDataAction==1 && IsUnrepresentedData(runNumber,lbn.second.second)) continue; //skip over bad lumi
+ lumisum += lbn.second.first;
+ if(lumisum >= lumi) return lbn.first;
}
-
- //the total lumi is stored in the run=0 slot of this map
- double lumi = myMap2->second[0] * m_random3->Rndm();
- double lumisum = 0.;
- std::map<UInt_t, Double_t>::const_iterator itend = myMap2->second.end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = myMap2->second.begin(); it2 != itend; ++it2)
- {
- //skip case of it2->first=0 (which is the total lumi)
- if(it2->first==0) continue;
- lumisum += it2->second;
- if (lumisum >= lumi)
- return it2->first;
- }
- Error("GetRandomRunNumber","overran integrated luminosity for mcRunNumber=%d, mubin=%d",mcRunNumber,muBin);
- throw std::runtime_error("Throwing 46: overran integrated luminosity for mcRunNumber");
- return 0;
+ Error("GetRandomLumiBlockNumber","overran integrated luminosity for RunNumber=%d",runNumber);
+ throw std::runtime_error("Throwing 46: overran integrated luminosity for runNumber");
+ return 0;
}
-Int_t Root::TPileupReweighting::GetRandomPeriodNumber(Int_t mcRunNumber) {
- if(!m_isInitialized) { Error("GetRandomPeriodNumber","Tool is not initialized"); return 0;}
+//only considers periods assigned directly!
+Int_t CP::TPileupReweighting::GetRandomPeriodNumber(Int_t periodNumber) {
+ if(!m_isInitialized) { Info("GetRandomPeriodNumber","Initializing the subtool.."); Initialize(); }
if(m_countingMode) { return 0; } //do nothing when in counting mode
- //check if the runNumber has been reassigned
- if(m_mcRemappings.find(mcRunNumber)!=m_mcRemappings.end()) mcRunNumber=m_mcRemappings[mcRunNumber];
-
- if(m_periods.find(mcRunNumber)==m_periods.end()) {
- Error("GetRandomPeriodNumber","Unknown MC RunNumber: %d",mcRunNumber);
- throw std::runtime_error("Throwing 46: Unknown MC RunNumber");
- }
- //if this runnumber has only one period, return itself
- if(m_periods[mcRunNumber].size()==1) {
- return mcRunNumber;
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetRandomPeriodNumber","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- const std::map<Int_t, Double_t>::const_iterator it = periodTotals["pileup"][-1].find(mcRunNumber);
- if(it==periodTotals["pileup"][-1].end()) {
- if(m_SetWarnings) Warning("GetRandomPeriodNumber","%d has no associated lumi. Returning 0",mcRunNumber);
- return 0; /*throw 33;*/
+ //if no subperiods, we just return ourselves
+ if(p->subPeriods.size()==0) return p->id;
+
+ //need the total lumi of this period
+ double lumi = p->sumOfWeights[-1] * m_random3->Rndm();
+
+ //loop over subPeriods ..
+ double lumiSum=0;
+ for(auto subp : p->subPeriods) {
+ lumiSum += subp->sumOfWeights[-1];
+ if(lumiSum >= lumi) return subp->id;
}
- Double_t lumi = it->second * m_random3->Rndm();
- Double_t lumisum = 0.;
- for(std::map<Int_t, std::pair<UInt_t,UInt_t> >::const_iterator it2 = m_periods[mcRunNumber].begin(); it2 != m_periods[mcRunNumber].end(); ++it2) {
- lumisum += periodTotals["pileup"][-1][it2->first];
- if (lumisum >= lumi)
- return it2->first;
- }
- Error("GetRandomPeriodNumber","overran integrated luminosity for mcRunNumber=%d",mcRunNumber);
- throw std::runtime_error("Throwing 46: overran integrated luminosity for mcRunNumber");
- return 0;
+
+ Error("GetRandomPeriodNumber","overran integrated luminosity for periodNumber=%d",periodNumber);
+ throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomPeriodNumber");
+ return 0;
+
}
-UInt_t Root::TPileupReweighting::GetRandomLumiBlockNumber(UInt_t runNumber) {
- if(!m_isInitialized) { Error("GetRandomRunNumber","Tool is not initialized"); return 0;}
- if(m_countingMode) { return 0; } //do nothing when in counting mode
+Bool_t CP::TPileupReweighting::IsUnrepresentedData(Int_t runNumber, Float_t x, Float_t y) {
+ int bin = m_emptyHistogram->FindFixBin(x,y);
+ return ( m_runs[runNumber].badBins[bin] );
+}
- double lumi = GetIntegratedLumi(runNumber,runNumber) * m_random3->Rndm();
- double lumisum = 0.;
- std::map<UInt_t, std::map<UInt_t, Double_t> >::iterator it = m_lumiByRunAndLbn.find(runNumber);
- if(it==m_lumiByRunAndLbn.end()) {
- Error("GetRandomLumiBlockNumber","Unknown run number: %d",runNumber);
- throw std::runtime_error("Throwing 46.1: Unknown Run number");
+//this method builds a file get can be friended to a TTree with a prw-hash branch
+Bool_t CP::TPileupReweighting::MakeWeightTree(TString channelNumbers, TString outFile, TString hashBranch, TString weightBranch) {
+ if(!m_isInitialized) { Info("MakeWeightTree","Initializing the subtool.."); Initialize(); }
+ TH1* hist = m_emptyHistogram;
+ if(!hist) {
+ Error("MakeWeightTree","Tool not configured properly ... please report this!");
+ throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
}
- std::map<UInt_t, Double_t>::const_iterator itend = it->second.end();
- for(std::map<UInt_t, Double_t>::const_iterator it2 = it->second.begin(); it2 != itend; ++it2) {
- lumisum += it2->second;
- if (lumisum >= lumi) return it2->first;
+
+ //loop over given channels, and loop over all periods (except the global period) and get the pileup weight in each case for each bin
+ TFile f1(outFile,"RECREATE");
+
+ TTree* outTree = new TTree("prwTree","prwTree");
+ ULong64_t prwHash(0);Float_t weight(0.);
+ outTree->Branch(hashBranch,&prwHash);
+ outTree->Branch(weightBranch,&weight);
+
+ TObjArray *tx = channelNumbers.Tokenize(",");
+ for (Int_t i = 0; i < tx->GetEntries(); i++) {
+ int channelNumber = ((TObjString *)(tx->At(i)))->String().Atoi();
+ if(channelNumber>999999 || channelNumber<0) {
+ Error("MakeWeightTree","ChannelNumber can not be bigger than 999999 or less than 0 ... got %d",channelNumber);
+ f1.Close();
+ return 0;
+ }
+ for(auto& period : m_periods) {
+ if(period.first==-1) continue;
+ int periodNumber = period.first;
+ for(int i=1;i<=hist->GetNbinsX();i++) {
+ double x = hist->GetXaxis()->GetBinCenter(i);
+ for(int j=1;j<=hist->GetNbinsY();j++) {
+ double y = hist->GetYaxis()->GetBinCenter(j);
+ weight = GetCombinedWeight(periodNumber,channelNumber,x,y);
+ prwHash = GetPRWHash(periodNumber,channelNumber,x,y);
+ outTree->Fill();
+ }
+ }
+
+ }
}
- Error("GetRandomLumiBlockNumber","overran integrated luminosity for RunNumber=%d",runNumber);
- throw std::runtime_error("Throwing 46: overran integrated luminosity for runNumber");
- return 0;
+
+ outTree->BuildIndex(hashBranch.Data());
+ outTree->Write();
+ f1.Close();
+
+ Info("MakeWeightTree","Successfully wrote prwTree to %s",outFile.Data());
+
+ return true;
}
-Float_t Root::TPileupReweighting::getPileupWeight(Float_t mu, Int_t runNumber, Int_t channel ) {
- if(!m_isInitialized) { Error("getPileupWeight","Tool is not initialized"); return 0;}
- if(m_SetWarnings) Warning("getPileupWeight","Depricated. Please use GetCombinedWeight(runNumber,channelNumber,mu)");
- return GetCombinedWeight(runNumber,channel,mu);
+ULong64_t CP::TPileupReweighting::GetPRWHash(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y) {
+ if(!m_isInitialized) { Info("GetPRWHash","Initializing the subtool.."); Initialize(); }
+ TH1* hist = m_emptyHistogram;
+ if(!hist) {
+ Error("GetPRWHash","Tool not configured properly ... please report this!");
+ throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
+ }
+
+ ULong64_t out = hist->FindFixBin(x,y);
+ out += (unsigned long)(channelNumber)*10000000000;
+ out += (unsigned long)(periodNumber)*10000;
+
+ return out;
+
}
-Float_t Root::TPileupReweighting::GetCombinedWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y, Float_t /*z*/) {
- if(!m_isInitialized) { Error("GetCombinedWeight","Tool is not initialized"); return 0;}
+Float_t CP::TPileupReweighting::GetCombinedWeight(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y) {
+ if(!m_isInitialized) { Info("GetCombinedWeight","Initializing the subtool.."); Initialize(); }
if(m_countingMode) return 0.;
if(m_countingMode) {
Error("GetCombinedWeight","You cannot get a weight when in config file generating mode. Please use FillWeight method");
throw std::runtime_error("Throwing 2: You cannot get a weight when in config file generating mode. Please use FillWeight method");
}
//decide how many dimensions this weight has - use the emptyHistogram to tell...
- TH1* hist = m_emptyHistograms[weightName];
+ TH1* hist = m_emptyHistogram;
if(!hist) {
- Error("GetCombinedWeight","Weight not recognised: %s",weightName.Data());
+ Error("GetCombinedWeight","Tool not configured properly ... please report this!");
+ throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
}
- //get the PeriodWeight for the associated runNumber - NOTE: Delayed implementing ... see Recipe H... check dependence... 09/11/2012 - looks like H was wrong to apply IntegratedLumiFraction and old CombinedWeight... since old CombinedWeight used the periodNumber's periodWeight, which appears to include the fraction already
- Int_t mcRunNumber = periodNumber;
- std::map<Int_t,Int_t>::iterator it = m_periodToMCRun.find(periodNumber);
- if(it!=m_periodToMCRun.end()) mcRunNumber = it->second;
- Float_t out = GetPeriodWeight(weightName,mcRunNumber/*periodNumber*/,channelNumber)*GetPrimaryWeight(weightName,periodNumber,channelNumber,x);
- if(hist->GetDimension()>1) out *= GetSecondaryWeight(weightName,periodNumber,channelNumber,x,y);
- //if(hist->GetDimension()>2) out *= GetTertiaryWeight(weightName,periodNumber,channelNumber,x,y,z);
-
+ Float_t out = GetPeriodWeight(periodNumber,channelNumber)*GetPrimaryWeight(periodNumber,channelNumber,x);
+ if(hist->GetDimension()>1) out *= GetSecondaryWeight(periodNumber,channelNumber,x,y);
+
return out;
}
-Float_t Root::TPileupReweighting::GetPeriodWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber) {
- if(!m_isInitialized) { Error("GetPeriodWeight","Tool is not initialized"); return 0;}
+Float_t CP::TPileupReweighting::GetPeriodWeight(Int_t periodNumber, Int_t channelNumber) {
+ //= L_A/L / N_A/N
+ if(!m_isInitialized) { Info("GetPeriodWeight","Initializing the subtool.."); Initialize(); }
if(m_countingMode) return 0.;
- //check if the runNumber has been reassigned
- //only do this for mc, never for data
- Int_t inPeriodNumber = periodNumber; //hold for getting the default channel if necessary
- if(channelNumber>=0) {
- if(m_mcRemappings.find(periodNumber)!=m_mcRemappings.end()) periodNumber=m_mcRemappings[periodNumber];
- }
- if(periodTotals[weightName][-1].find(periodNumber)==periodTotals[weightName][-1].end()) {
- Error("GetPeriodWeight","Unrecognised periodNumber/mcRunNumber: %d",periodNumber);
- throw std::runtime_error("Throwing 1: Unrecognised periodNumber/mcRunNumber");
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetPeriodWeight","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
+ if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) channelNumber = p->defaultChannel;
- Double_t data = periodTotals[weightName][-1][periodNumber];
- Double_t totalData = globalTotals[weightName][-1];
+ double n_a = p->sumOfWeights[channelNumber];
+ double n = m_periods[-1]->sumOfWeights[channelNumber];
- if(totalData==0.) {
- Error("GetPeriodWeight","No data for weight %s",weightName.Data());
- throw std::runtime_error("Throwing 2: No data for weight");
- }
+ double l_a = p->sumOfWeights[-1];
+ double l = m_periods[-1]->sumOfWeights[-1];
- if(periodTotals[weightName].find(channelNumber)==periodTotals[weightName].end()) {
- //use the default
- if(periodTotals[weightName].find(GetDefaultChannel(inPeriodNumber))==periodTotals[weightName].end()) {
- Error("GetPeriodWeight","Unrecognised channel %d and no default found: %d",channelNumber,GetDefaultChannel(inPeriodNumber));
- throw std::runtime_error("Throwing 3: Unrecognised channel");
- }
- channelNumber=GetDefaultChannel(inPeriodNumber);
- }
- if(periodTotals[weightName][channelNumber].find(periodNumber)==periodTotals[weightName][channelNumber].end()) {
- Error("GetPeriodWeight","Unrecognised periodNumber/mcRunNumber %d in channelNumber %d",periodNumber,channelNumber);
- throw std::runtime_error("Throwing 1: Unrecognised periodNumber/mcRunNumber");
- }
+ return (l_a/l) / (n_a/n);
+}
- Double_t mc = periodTotals[weightName][channelNumber][periodNumber];
- Double_t totalMC = globalTotals[weightName][channelNumber];
+Float_t CP::TPileupReweighting::GetPrimaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x) {
+ //= L_i/L_A / N_i/N_A .. primaryHists have already been normalized
+ if(!m_isInitialized) { Info("GetPrimaryWeight","Initializing the subtool.."); Initialize(); }
+ if(m_countingMode) return 0.;
- if(totalMC==0.) {
- if(m_SetWarnings) Warning("GetPeriodWeight","No mc for weight '%s' in channelNumber %d",weightName.Data(),channelNumber);
- return 0.;
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetPrimaryWeight","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- if(mc==0.) {
- if(m_SetWarnings) Warning("GetPeriodWeight","No mc for weight '%s' in channelNumber %d and periodNumber/runNumber %d",weightName.Data(),channelNumber,periodNumber);
- return 0.;
+ int oChanNumber = channelNumber;
+ if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) channelNumber = p->defaultChannel;
+
+ if(!p->primaryHists[channelNumber]) {
+ Error("GetPrimaryWeight","Unrecognised channelNumber: %d",oChanNumber);
+ throw std::runtime_error("Throwing 2: Unrecognised channelNumber");
}
- return (data/totalData) / (mc/totalMC) ;
+ int bin = p->primaryHists[channelNumber]->FindFixBin(x);
-}
+ double n = p->primaryHists[channelNumber]->GetBinContent(bin);
-Float_t Root::TPileupReweighting::GetPrimaryWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber, Float_t x) {
- if(!m_isInitialized) { Error("GetPrimaryWeight","Tool is not initialized"); return 0;}
- if(m_countingMode) return 0.;
- //check if the runNumber has been reassigned
- //only do this for mc, never for data
- Int_t inPeriodNumber = periodNumber;
- if(channelNumber>=0) {
- if(m_mcRemappings.find(periodNumber)!=m_mcRemappings.end()) periodNumber=m_mcRemappings[periodNumber];
+ if(!p->primaryHists[-1]) {
+ Error("GetPrimaryWeight","No data loaded for period %d. Did you forget to load a lumicalc file or data config file?",periodNumber);
+ throw std::runtime_error("Throwing 3: No data loaded. Did you forget to load a lumicalc file or data config file?");
}
- //determine bin number
- TH1D* data = primaryDistributions[weightName][-1][periodNumber];
- if(!data) {
- Error("GetPrimaryWeight","Unrecognised periodNumber/mcRunNumber: %d",periodNumber);
- throw std::runtime_error("Throwing 4: Unrecognised periodNumber/mcRunNumber");
- }
- if(primaryDistributions[weightName].find(channelNumber)==primaryDistributions[weightName].end()) {
- //try switching to default channel
- if(primaryDistributions[weightName].find(GetDefaultChannel(inPeriodNumber))==primaryDistributions[weightName].end()) {
- Error("GetPrimaryWeight","Unrecognised channelNumber %d, and no distribution found for defaultChannel %d",channelNumber,GetDefaultChannel(inPeriodNumber));
- throw std::runtime_error("Throwing 5: Unrecognised channelNumber");
- }
- channelNumber=GetDefaultChannel(inPeriodNumber);
- }
- TH1D* mc = primaryDistributions[weightName][channelNumber][periodNumber];
- if(!mc) {
- Error("GetPrimaryWeight","Unable to find distribution for [%s,%d,%d]",weightName.Data(),channelNumber,periodNumber);
- throw std::runtime_error("Throwing 6: Unable to find distribution");
- }
+ double l = p->primaryHists[-1]->GetBinContent(bin);
- //MC10b correction
- if ( fabs(x - 100.0) < 0.00001 && mc->GetXaxis()->GetXmax() < 99.0 ) x = 0.0;
+ return l/n;
+}
- Int_t bin = data->GetXaxis()->FindFixBin(x);
+Float_t CP::TPileupReweighting::GetSecondaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x,Float_t y) {
+ //= L_j/L_i / N_j/N_i .. secondary hists have already been normalized
+ if(!m_isInitialized) { Info("GetSecondaryWeight","Initializing the subtool.."); Initialize(); }
+ if(m_countingMode) return 0.;
- if(mc->GetBinContent(bin)==0.) {
- if(m_SetWarnings) Warning("GetPrimaryWeight","No mc for weight '%s' in channelNumber=%d, periodNumber/runNumber=%d,x=%f",weightName.Data(),channelNumber,periodNumber,x);
- return 0.;
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetSecondaryWeight","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
+ if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) channelNumber = p->defaultChannel;
+ int bin = p->secondaryHists[channelNumber]->FindFixBin(x,y);
+ double n = p->secondaryHists[channelNumber]->GetBinContent(bin);
+ double l = p->secondaryHists[-1]->GetBinContent(bin);
+
+ return l/n;
+}
- return (data->GetBinContent(bin)/mc->GetBinContent(bin));
+Double_t CP::TPileupReweighting::GetDataWeight(Int_t runNumber, TString trigger) {
+ //special mu-independent version of GetDataWeight. Will just use the global luminosity
+ m_doGlobalDataWeight=true;
+ double out = GetDataWeight(runNumber, trigger, 0);
+ m_doGlobalDataWeight=false;
+ return out;
}
-Float_t Root::TPileupReweighting::GetSecondaryWeight(const TString weightName, Int_t periodNumber, Int_t channelNumber,Float_t x, Float_t y) {
- if(!m_isInitialized) { Error("GetSecondaryWeight","Tool is not initialized"); return 0;}
+Double_t CP::TPileupReweighting::GetDataWeight(Int_t runNumber, TString trigger, Double_t x) {
+
+ if(!m_isInitialized) { Info("GetDataWeight","Initializing the subtool.."); Initialize(); }
if(m_countingMode) return 0.;
- //check if the runNumber has been reassigned
- //only do this for mc, never for data
- Int_t inPeriodNumber = periodNumber;
- if(channelNumber>=0) {
- if(m_mcRemappings.find(periodNumber)!=m_mcRemappings.end()) periodNumber=m_mcRemappings[periodNumber];
+ if(m_countingMode) {
+ Error("GetDataWeight","You cannot get a weight when in config file generating mode. Please use FillWeight method");
+ throw std::runtime_error("GetDataWeight: You cannot get a weight when in config file generating mode. Please use FillWeight method");
}
+ //determine which period this run number is in
+ Int_t periodNumber = GetFirstFoundPeriodNumber(runNumber);
- //determine bin number
- TH2D* data = secondaryDistributions[weightName][-1][periodNumber];
- if(!data) {
- Error("GetSecondaryWeight","Unrecognised periodNumber/mcRunNumber: %d",periodNumber);
- throw std::runtime_error("Throwing 4: Unrecognised periodNumber/mcRunNumber");
+ Period* p = m_periods[periodNumber];
+ if(!p) {
+ Error("GetDataWeight","Unrecognised periodNumber: %d",periodNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- if(secondaryDistributions[weightName].find(channelNumber)==secondaryDistributions[weightName].end()) {
- //try switching to default channel
- if(secondaryDistributions[weightName].find(GetDefaultChannel(inPeriodNumber))==secondaryDistributions[weightName].end()) {
- Error("GetSecondaryWeight","Unrecognised channelNumber %d, and no distribution found for defaultChannel %d",channelNumber,GetDefaultChannel(inPeriodNumber));
- throw std::runtime_error("Throwing 5: Unrecognised channelNumber");
- }
- channelNumber=GetDefaultChannel(inPeriodNumber);
+
+ //see if we already made this trigger hist
+
+
+ if(p->triggerHists.find(trigger)==p->triggerHists.end()) {
+ //try to do a subtrigger calculation we need to reopen all the lumicalc files and do the fiddily prescaled luminosity calculation
+ //will throw errors if can't
+ CalculatePrescaledLuminosityHistograms(trigger);
}
- TH2D* mc = secondaryDistributions[weightName][channelNumber][periodNumber];
- if(!mc) {
- Error("GetSecondaryWeight","Unable to find distribution for [%s,%d,%d]",weightName.Data(),channelNumber,periodNumber);
- throw std::runtime_error("Throwing 6: Unable to find distribution");
+ TH1D* denomHist = p->triggerHists[trigger];
+ if(denomHist==0) {
+ Error("GetDataWeight","Could not find trigger %s in period %d",trigger.Data(),p->id);
+ throw std::runtime_error("GetDataWeight: Could not find trigger");
}
- Int_t nx = data->GetXaxis()->GetNbins()+2;
- Int_t binx = data->GetXaxis()->FindFixBin(x);
- Int_t biny = data->GetYaxis()->FindFixBin(y);
- Int_t bin = binx + nx*biny;
-
- if(mc->GetBinContent(bin)==0.) {
- if(m_SetWarnings) Warning("GetSecondaryWeight","No mc for weight '%s' in channelNumber %d and periodNumber/runNumber %d",weightName.Data(),channelNumber,periodNumber);
- return 0.;
+ TH1D* numerHist = p->triggerHists["None"];
+ if(numerHist==0) {
+ Error("GetDataWeight","Could not find unprescaled trigger in period %d",p->id);
+ throw std::runtime_error("GetDataWeight: Could not find unprescaled trigger. Please AddLumiCalc with a 'None' trigger");
}
+
+ if(m_doGlobalDataWeight) return numerHist->Integral(0,numerHist->GetNbinsX()+1)/denomHist->Integral(0,denomHist->GetNbinsX()+1);
- return (data->GetBinContent(bin)/mc->GetBinContent(bin));
-
-}
+ Int_t bin=numerHist->FindFixBin(x);
+ return numerHist->GetBinContent(bin)/denomHist->GetBinContent(bin);
-Int_t Root::TPileupReweighting::Fill(Int_t runNumber,Int_t channelNumber,Float_t w, Float_t x, Float_t y, Float_t z) {
- return Fill("pileup",runNumber,channelNumber,w,x,y,z);
}
//fills the appropriate inputHistograms
-Int_t Root::TPileupReweighting::Fill(const TString weightName,Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y, Float_t z) {
- if(!m_isInitialized) { Error("Fill","Tool is not initialized!"); return 0;}
+Int_t CP::TPileupReweighting::Fill(Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y) {
+ if(!m_isInitialized) { Info("Fill","Initializing the subtool.."); Initialize(); }
//should only be given genuine mcRunNumbers if mc (channel>=0). We don't fill periodNumber distributions
- if(channelNumber>=0) {
- //if(m_mcRemappings.find(runNumber)!=m_mcRemappings.end()) runNumber=m_mcRemappings[runNumber]; //- no remappings in fill mode!
- //also check that channel histograms exist for all defined periods. If not, define them
- if(m_inputHistograms[weightName].find(channelNumber)==m_inputHistograms[weightName].end()) {
- //loop over mcruns, get empty histogram for each
- for(std::map<Int_t,std::map<Int_t, std::pair<UInt_t,UInt_t> > >::iterator it=m_periods.begin();it!=m_periods.end();++it) {
- m_inputHistograms[weightName][channelNumber][it->first] = CloneEmptyHistogram(weightName,it->first,channelNumber);
- }
- //must create histograms for the merged run numbers as well, as config files are made without merging the numbers
- //only do it if the remapping has some actual period assignments though!
- for(std::map<Int_t,Int_t>::iterator it=m_mcRemappings.begin();it!=m_mcRemappings.end();++it) {
- if(m_periods.find(it->second) == m_periods.end()) continue; ///added this line in july 2014... need to check doesnt break MC12ab
- m_inputHistograms[weightName][channelNumber][it->first] = CloneEmptyHistogram(weightName,it->first,channelNumber);
- }
+ TH1* hist = 0;
+
+ //get the period if mc, get the data run if data
+ if(channelNumber>=0) {
+ if(m_periods.find(runNumber)==m_periods.end()) {
+ Error("Fill","Unrecognised runNumber: %d. Check your period configuration (AddPeriod or UsePeriodConfig)",runNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
- } else {
- //create histograms for data too, only for this particular run
- if(m_inputHistograms[weightName][channelNumber].find(runNumber)==m_inputHistograms[weightName][channelNumber].end()) {
- m_inputHistograms[weightName][channelNumber][runNumber] = CloneEmptyHistogram(weightName,runNumber,channelNumber);
+ Period* p = m_periods[runNumber];
+ if(!p) {
+ Error("Fill","Unrecognised runNumber: %d. Check your period configuration (AddPeriod or UsePeriodConfig) ... but should never have got here so please report this!",runNumber);
+ throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
}
+ if(p->inputHists.find(channelNumber)==p->inputHists.end()) {
+ //need to create my period histogram
+ p->inputHists[channelNumber] = CloneEmptyHistogram(runNumber,channelNumber);
+ }
+ hist = p->inputHists[channelNumber];
+ } else {
+ Run& r = m_runs[runNumber];
+ if(r.inputHists["None"]==0) r.inputHists["None"]=CloneEmptyHistogram(runNumber,channelNumber);
+ hist = r.inputHists["None"];
}
- TH1* hist = m_inputHistograms[weightName][channelNumber][runNumber];
-
if(!hist) {
- Error("Fill","Unknown [weight,run,channel] = [%s,%d,%d]",weightName.Data(),runNumber,channelNumber);
- throw std::runtime_error("Throwing 45: Fill:: Unknown [weight,run,channel] ");
+ Error("Fill","Unknown [run,channel] = [%d,%d]",runNumber,channelNumber);
+ throw std::runtime_error("Throwing 45: Fill:: Unknown [run,channel] ");
}
if(hist->GetDimension()==1) {
return hist->Fill(x,w);
} else if(hist->GetDimension()==2) {
- return (dynamic_cast<TH2*>(hist))->Fill(x,y,w);
- } else if(hist->GetDimension()==3) {
- return (dynamic_cast<TH3*>(hist))->Fill(x,y,z,w);
- }
+ return (static_cast<TH2*>(hist))->Fill(x,y,w);
+ }
return -1;
}
+Int_t CP::TPileupReweighting::WriteToFile(const TString fname) {
+
+ if(!m_countingMode) {Warning("WriteToFile","Not in counting mode, so no file will be written");return 0;}
+
+
+ //build a TTree with the correct mc structure, and dump the mc histogram info in to it
+ //also build aggregate histograms across all channels - shouldn't be used as input histograms unless for a single channel
+ TString filename = fname;
+ filename += (filename=="") ? TString(this->GetName()) + ".prw.root" : "";
+
+ //loop over the weights. "pileup" gets it own MCPileupReweighting ttree. ... since we got rid of weightNames in this version, everything goes in MCPileupReweighting now!
+ //data goes in to DataCustomReweighting. other goes in to MCCustomReweighting
+ TFile* outFile = TFile::Open(filename,"RECREATE");
+ Int_t r = WriteToFile(outFile);
+ outFile->Close();
+ return r;
+}
+
+
+Int_t CP::TPileupReweighting::WriteToFile(TFile* outFile) {
+ if(!m_countingMode) {Warning("WriteToFile","Not in counting mode, so no file will be written");return 0;}
+
+ TDirectory* origDir = gDirectory;
+
+ outFile->cd();
+
+ TTree *outTreeMC=0;
+ TTree *outTreeData=0;
+ Int_t channel = 0;UInt_t runNumber = 0;
+ std::vector<UInt_t>* pStarts = 0;std::vector<UInt_t>* pEnds = 0;
+ Char_t histName[150];
+
+
+ //loop over periods ... periods only get entry in table if they have an input histogram
+ for(auto period : m_periods) {
+ if(period.first != period.second->id) continue; //skips redirects
+ runNumber = period.first;
+ pStarts = new std::vector<UInt_t>;pEnds = new std::vector<UInt_t>;
+ if(period.second->subPeriods.size()==0) {
+ pStarts->push_back(period.second->start); pEnds->push_back(period.second->end);
+ }
+ else {
+ for(auto subp : period.second->subPeriods) {
+ pStarts->push_back(subp->start); pEnds->push_back(subp->end);
+ }
+ }
+ for(auto inHist : period.second->inputHists) {
+ channel = inHist.first;
+ TH1* hist = inHist.second;
+ strcpy(histName,hist->GetName());
+ hist->Write();
+
+ if(!outTreeMC) {
+ outTreeMC = new TTree("MCPileupReweighting","MCPileupReweighting");
+ outTreeMC->Branch("Channel",&channel);
+ outTreeMC->Branch("RunNumber",&runNumber);
+ outTreeMC->Branch("PeriodStarts",&pStarts);
+ outTreeMC->Branch("PeriodEnds",&pEnds);
+ outTreeMC->Branch("HistName",&histName,"HistName/C");
+ }
+ outTreeMC->Fill();
+ }
+ delete pStarts; delete pEnds;
+ }
+
+ //loop over data
+ for(auto& run : m_runs) {
+ runNumber = run.first;
+ if(run.second.inputHists.find("None")==run.second.inputHists.end()) continue;
+
+ TH1* hist = run.second.inputHists["None"];
+ strcpy(histName,hist->GetName());
+ hist->Write();
+ if(!outTreeData) {
+ outTreeData = new TTree("DataPileupReweighting","DataPileupReweighting");
+ outTreeData->Branch("RunNumber",&runNumber);
+ outTreeData->Branch("HistName",&histName,"HistName/C");
+ }
+ outTreeData->Fill();
+ }
+
+
+ //write the non-zero ttrees
+ if(outTreeMC) {outTreeMC->Write();delete outTreeMC;outTreeMC=0;}
+ if(outTreeData) {outTreeData->Write();delete outTreeData;outTreeData=0;}
+
+
+ Info("WriteToFile", "Successfully generated config file: %s",outFile->GetName());
+ Info("WriteToFile", "Happy Reweighting :-)");
+
+ gDirectory = origDir;
+
+ return 0;
+}
+
+
//=============================================================================
// Method to calculate the ratio histogram
//=============================================================================
-void Root::TPileupReweighting::normalizeHistogram(TH1* hist){
+void CP::TPileupReweighting::normalizeHistogram(TH1* hist){
// normalize the data histogram based on which sort of histogram it is
if(hist){
if(hist->InheritsFrom("TH3")) {
@@ -2064,7 +1645,7 @@ void Root::TPileupReweighting::normalizeHistogram(TH1* hist){
//This allows PROOF to merge the generated histograms before the WriteToFile or WriteCustomToFile calls
-Int_t Root::TPileupReweighting::Merge(TCollection *coll) {
+Int_t CP::TPileupReweighting::Merge(TCollection *coll) {
if(!coll) return 0;
if(coll->IsEmpty()) return 0;
@@ -2074,27 +1655,32 @@ Int_t Root::TPileupReweighting::Merge(TCollection *coll) {
while( ( obj = next() ) ) {
// Check that the element is of the right type:
- Root::TPileupReweighting* vobj = dynamic_cast< Root::TPileupReweighting* >( obj );
+ CP::TPileupReweighting* vobj = dynamic_cast< CP::TPileupReweighting* >( obj );
if( ! vobj ) {
Error( "Merge", "Unknown object type encountered: %s",obj->ClassName() );
return 0;
}
- //merge the inputHistograms
- for(std::map<TString, std::map<Int_t,std::map<Int_t, TH1*> > >::iterator weights=vobj->GetInputHistograms().begin();weights!=vobj->GetInputHistograms().end();++weights) {
- TString weightName = weights->first;
- for(std::map<Int_t,std::map<Int_t, TH1*> >::iterator channels=weights->second.begin();channels!=weights->second.end();++channels) {
- Int_t channel = channels->first;
- for(std::map<Int_t, TH1*>::iterator runs=channels->second.begin();runs!=channels->second.end();++runs) {
- Int_t run=runs->first;
- if(runs->second) {
- if(m_inputHistograms[weightName][channel][run]) {
- m_inputHistograms[weightName][channel][run]->Add(runs->second);
- } else {
- m_inputHistograms[weightName][channel][run] = dynamic_cast<TH1*>(runs->second->Clone(runs->second->GetName()));
- m_inputHistograms[weightName][channel][run]->SetDirectory(0);
- }
- }
+ //merge the inputHistograms ... all the periods should be identical
+ for(auto period : vobj->m_periods) {
+ if(period.first != period.second->id) continue;
+ for(auto& iHist : period.second->inputHists) {
+ if(GetInputHistogram(period.first,iHist.first)==0) {
+ m_periods[period.first]->inputHists[iHist.first] = dynamic_cast<TH1*>(iHist.second->Clone(iHist.second->GetName()));
+ m_periods[period.first]->inputHists[iHist.first]->SetDirectory(0);
+ } else {
+ GetInputHistogram(period.first,iHist.first)->Add(iHist.second);
+ }
+ }
+ }
+ //also must remember to merge the runs too (where the data is held)
+ for(auto& run : vobj->m_runs) {
+ for(auto& iHist : run.second.inputHists) {
+ if(m_runs[run.first].inputHists[iHist.first]==0) {
+ m_runs[run.first].inputHists[iHist.first] = dynamic_cast<TH1*>(iHist.second->Clone(iHist.second->GetName()));
+ m_runs[run.first].inputHists[iHist.first]->SetDirectory(0);
+ } else {
+ m_runs[run.first].inputHists[iHist.first]->Add(iHist.second);
}
}
}
@@ -2104,7 +1690,7 @@ Int_t Root::TPileupReweighting::Merge(TCollection *coll) {
}
-void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TString trigger) {
+void CP::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TString trigger) {
//check if this is a composite trigger (user has provided OR of triggers that are passed before prescale)
TString triggerCopy = trigger; triggerCopy.ReplaceAll(" ",""); triggerCopy.ReplaceAll("&&","&");triggerCopy.ReplaceAll("||","|");
@@ -2116,25 +1702,6 @@ void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TStr
std::vector<TString> subTriggers;
-/*
- int fStart=0;int nextSep=trigger.Index(',');
- while(nextSep>=0) {
- TString subTrigger = trigger(fStart,nextSep-fStart);
- if(m_triggerPrimaryDistributions.find(subTrigger)==m_triggerPrimaryDistributions.end()) {
- Error("GetDataWeight","Could not find subTrigger %s in composite trigger set: %s",subTrigger.Data(),trigger.Data());
- throw std::runtime_error("Could not find sub-trigger");
- }
- subTriggers.push_back(subTrigger);
- fStart = nextSep+1;
- nextSep = trigger.Index(',',fStart);
- }
- subTriggers.push_back(trigger(fStart,trigger.Length()-fStart)); //don't forget the last one
- if(subTriggers.size()==0) {
- Error("GetDataWeight","Could not find trigger %s",trigger.Data());
- throw std::runtime_error("Could not find trigger");
- }
-*/
-
t->getTriggers(subTriggers); //fills the vector
//1. Open all the lumicalc files
@@ -2180,6 +1747,7 @@ void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TStr
long nEntries = tmp->GetEntries();
for(long i=0;i<nEntries;i++) {
b_runNbr->GetEntry(i);b_L1Presc->GetEntry(i);b_L2Presc->GetEntry(i);b_L3Presc->GetEntry(i);b_lbn->GetEntry(i);
+ runNbr += m_lumicalcRunNumberOffset;
//save the prescale by run and lbn
//if(runNbr==215643) Info("...","prescale in [%d,%d] = %f %f %f", runNbr,lbn,ps1,ps2,ps3);
if(ps1>0&&ps2>0&&ps3>0) prescaleByRunAndLbn[*it][runNbr][lbn] = ps1*ps2*ps3;
@@ -2220,6 +1788,7 @@ void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TStr
b_runNbr->GetEntry(i);b_intLumi->GetEntry(i);b_mu->GetEntry(i);
b_lbn->GetEntry(i);
+ runNbr += m_lumicalcRunNumberOffset;
//for each subtrigger, lookup prescale, and calculate pFactor
/*
double pFactor(1);
@@ -2235,12 +1804,23 @@ void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TStr
//Info("...","prescale in [%d,%d] = %f", runNbr,lbn,1./pFactor);
- Int_t periodNumber = GetFirstFoundPeriodNumber(runNbr);
- if(m_triggerPrimaryDistributions[trigger].find(periodNumber) == m_triggerPrimaryDistributions[trigger].end()) {
- m_triggerPrimaryDistributions[trigger][periodNumber] = dynamic_cast<TH1D*>(CloneEmptyHistogram("pileup",periodNumber,-1));
- Info("CalculatePrescaledLuminosityHistograms","Created Data Weight Histogram for [%s,%d]",trigger.Data(),periodNumber);
+ int bin = m_emptyHistogram->FindFixBin(mu*m_dataScaleFactorX);
+
+ //add into all periods that contain this runNbr
+ for(auto p : m_periods) {
+ if(p.first != p.second->id) continue; //skips remappings
+ if(!p.second->contains(runNbr)) continue;
+ std::map<TString,TH1D*>& triggerHists = p.second->triggerHists;
+ if(triggerHists.find(trigger) == triggerHists.end()) {
+ triggerHists[trigger] = dynamic_cast<TH1D*>(CloneEmptyHistogram(p.first,-1));
+ if(m_debugging) Info("CalculatePrescaledLuminosityHistograms","Created Data Weight Histogram for [%s,%d]",trigger.Data(),p.first);
+ }
+ //check if we were about to fill a bad bin ... if we are, we either skipp the fill (unrep action=1) or redirect (unrep action=3)
+ if( (m_unrepresentedDataAction==1) && p.second->inputBinRedirect[bin]!=bin) { } //do nothing
+ else if( (m_unrepresentedDataAction==3) ) triggerHists[trigger]->Fill(triggerHists[trigger]->GetBinCenter(p.second->inputBinRedirect[bin]), intLumi*pFactor);
+ else triggerHists[trigger]->Fill(mu*m_dataScaleFactorX,intLumi*pFactor);
+
}
- m_triggerPrimaryDistributions[trigger][periodNumber]->Fill(mu*m_dataScaleFactorX,intLumi*pFactor);
}
}
}
@@ -2256,8 +1836,8 @@ void Root::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TStr
}
-Root::TPileupReweighting::CompositeTrigger* Root::TPileupReweighting::makeTrigger(TString s) {
- Info("makeTrigger","Doing %s",s.Data());
+CP::TPileupReweighting::CompositeTrigger* CP::TPileupReweighting::makeTrigger(TString s) {
+ if(m_debugging) Info("makeTrigger","Doing %s",s.Data());
//find the first operand
TString oper1; CompositeTrigger* cOper1 = 0;
if(s.BeginsWith("(")) {
@@ -2283,8 +1863,9 @@ Root::TPileupReweighting::CompositeTrigger* Root::TPileupReweighting::makeTrigge
i++;
}
+
if(op==0) {
- if(m_triggerPrimaryDistributions.find(s)==m_triggerPrimaryDistributions.end()) {
+ if(m_lumicalcFiles.find(s)==m_lumicalcFiles.end()) {
Error("GetDataWeight","Could not find subTrigger %s",s.Data());
return 0;
}
@@ -2314,12 +1895,12 @@ Root::TPileupReweighting::CompositeTrigger* Root::TPileupReweighting::makeTrigge
if(j==s.Length()+1 && op==2) { //no more OR found, set oper2 to remainder
if(cOper1) out->trig1=cOper1;
else {oper1=s(0,i); out->trig1=makeTrigger(oper1); }
- TString oper2 = s(i+1,s.Length()); Info("makeTrigger","Found & %s %s",oper1.Data(),oper2.Data());
+ TString oper2 = s(i+1,s.Length()); if(m_debugging) Info("makeTrigger","Found & %s %s",oper1.Data(),oper2.Data());
out->trig2 = makeTrigger(oper2); if(out->trig2==0) { delete out; return 0; }
return out;
} else if(op==1) { //found an OR, set oper1 to everything up to this
- oper1 = s(0,j); TString oper2 = s(j+1,s.Length()); Info("makeTrigger","Found & then | %s %s",oper1.Data(),oper2.Data());
+ oper1 = s(0,j); TString oper2 = s(j+1,s.Length()); if(m_debugging) Info("makeTrigger","Found & then | %s %s",oper1.Data(),oper2.Data());
if(cOper1) delete cOper1;
out->op=op; //updates to an OR
out->trig1=makeTrigger(oper1); if(out->trig1==0) { delete out; return 0; }
@@ -2331,28 +1912,5 @@ Root::TPileupReweighting::CompositeTrigger* Root::TPileupReweighting::makeTrigge
return 0;
}
-Double_t Root::TPileupReweighting::GetIntegratedLumi(TString trigger) {
- if(!m_isInitialized) {
- Error("GetIntegratedLumi", "Please initialize the tool before retrieving the total lumi");
- throw std::runtime_error("Throwing 1");
- }
- if(!m_lumiVectorIsLoaded) {
- Error("GetIntegratedLumi","No Lumicalc file loaded, so cannot get integrated lumi, returning 0");
- return 0;
- }
- if(trigger=="") return globalTotals["pileup"][-1]/1E6;
- else { //check the triggerPrimaryDistributions, if not present, build it
- if(m_triggerPrimaryDistributions.find(trigger)==m_triggerPrimaryDistributions.end()) {
- //try to do a subtrigger calculation we need to reopen all the lumicalc files and do the fiddily prescaled luminosity calculation
- //will throw errors if can't
- CalculatePrescaledLuminosityHistograms(trigger);
- }
- //loop over all periods and sum up lumi
- Double_t out =0;
- for(std::map<Int_t, TH1D*>::iterator it=m_triggerPrimaryDistributions[trigger].begin(); it!= m_triggerPrimaryDistributions[trigger].end(); ++it) {
- out += it->second->Integral(0,it->second->GetNbinsX()+1);
- }
- return out;
- }
-}
+
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/Makefile.RootCore b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/Makefile.RootCore
index 420f80bf502a8c365ddb8d8939e0ead9e206cc31..c88e9b711f76a356eb80ceaa823501ced71bc506 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/Makefile.RootCore
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/Makefile.RootCore
@@ -2,7 +2,8 @@ PACKAGE = PileupReweighting
PACKAGE_CXXFLAGS =
PACKAGE_LDFLAGS =
PACKAGE_PRELOAD = Hist
-PACKAGE_DEP = AsgTools xAODEventInfo PathResolver xAODCore
+PACKAGE_DEP = AsgTools ReweightUtils xAODEventInfo PathResolver xAODCore PATInterfaces
PACKAGE_REFLEX = 1
+PACKAGE_PEDANTIC = 1
include $(ROOTCOREDIR)/Makefile-common
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/requirements b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/requirements
index 074c81cfd491076054170ecc793966bdaadc0676..ee281aa266a6b82b48ab4479d99450ac9fcdc83e 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/requirements
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/cmt/requirements
@@ -9,7 +9,7 @@ use GaudiInterface GaudiInterface-* External
use AsgTools AsgTools-* Control/AthToolSupport
use xAODEventInfo xAODEventInfo-* Event/xAOD
-
+use PATInterfaces PATInterfaces-* PhysicsAnalysis/AnalysisCommon
private
use PathResolver PathResolver-* Tools
@@ -23,19 +23,31 @@ end_private
apply_tag ROOTBasicLibs
apply_tag ROOTMathLibs
+apply_tag ROOTCintexLibs
+
# Specify the required ROOT components for cmake (transparent to CMT)
apply_pattern cmake_add_command command="find_package(ROOT COMPONENTS MathCore RIO)"
## declare the directories CMT should know about
branches python share doc src Root
-apply_pattern have_root_headers root_headers="TPileupReweighting.h ../Root/LinkDef.h" headers_lib="PileupReweightingLib" extra_includes=" -I$(PathResolver_include_dirs) -I$(AsgTools_include_dirs) -DXAOD_STANDALONE"
-apply_pattern dual_use_library files=" *.cxx ../Root/*.cxx"
+
+library PileupReweightingLib "../Root/*.cxx"
+apply_pattern named_installed_library library=PileupReweightingLib
+
+macro extra_in "" AthAnalysisBase "-DXAOD_ANALYSIS"
+apply_pattern have_root_headers root_headers="TPileupReweighting.h ../Root/LinkDef.h" headers_lib="PileupReweightingLib" extra_includes="$(extra_in) "
+
+library PileupReweighting *.cxx -s=components *.cxx
+macro_append PileupReweighting_dependencies " PileupReweightingLib"
+apply_pattern component_library
+
+#apply_pattern dual_use_library files=" *.cxx ../Root/*.cxx"
private
use AtlasReflex AtlasReflex-* External -no_auto_imports
-apply_pattern lcgdict dict=PileupReweightingDict selectionfile=selection.xml headerfiles="../PileupReweighting/PileupReweightingDict.h"
+apply_pattern lcgdict dict=PileupReweighting selectionfile=selection.xml headerfiles="../PileupReweighting/PileupReweightingDict.h"
end_private
@@ -56,4 +68,6 @@ apply_pattern declare_calib files="../share/*.root"
-
+application testPRW testPRW.C
+#ensure the main library is compiled first
+macro_append testPRW_dependencies " PileupReweighting "
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/exampleJobOptions.py b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/exampleJobOptions.py
index f2c769d5ab0b7474c3f457584fc9ba91a51e543b..87f1b3e8be937b29d5879cf440328baa44ce1fe6 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/exampleJobOptions.py
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/exampleJobOptions.py
@@ -1,12 +1,27 @@
#by default (i.e. with no config files) the tool will be created in config file making mode
+from glob import glob
+
+theApp.EvtMax=10 #says how many events to run over. Set to -1 for all events
+
+import AthenaPoolCnvSvc.ReadAthenaPool #sets up reading of POOL files (e.g. xAODs)
+svcMgr.EventSelector.InputCollections=glob(vars().get("FILES","/afs/cern.ch/atlas/project/PAT/xAODs/r5591/mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591/AOD.01494882._111853.pool.root.1") ) #insert your list of input files here
+
+algseq = CfgMgr.AthSequencer("AthAlgSeq")
+
+
+#in your c++ code, create a ToolHandle<IPileupReweightingTool>
+#the ToolHandle constructor should be given "CP::PileupReweightingTool/myTool" as its string argument
ToolSvc += CfgMgr.CP__PileupReweightingTool("myTool",
- ConfigFiles=["PileupReweighting/mc12ab_defaults.prw.root"],
+ ConfigFiles=["my.prw.root"],
DataScaleFactor=1./1.09,
- UnrepresentedDataAction=2,
- LumiCalcFiles=["GRL_v61.None.lumicalc.root"])
+ UnrepresentedDataAction=2, DataScaleFactorUP=1./1.11, DataScaleFactorDOWN=1./1.07,
+ LumiCalcFiles=["/usera/will/testareas/ZZ8TeV/ZZAnalysis/share/GRL_v61.None.lumicalc.root"],OutputLevel=VERBOSE)
+
+algseq += CfgMgr.CP__PileupReweightingProvider(Tool=ToolSvc.myTool,OutputLevel=VERBOSE) #adds an instance of your alg to it
+
+
+include("AthAnalysisBaseComps/SuppressLogging.py") #Optional include to suppress as much athena output as possible
-#in your code, create a ToolHandle<IPileupReweightingTool>
-#the ToolHandle constructor should be given "CP::PileupReweightingTool/myTool" as its string argument
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/generatePRW_jobOptions.py b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/generatePRW_jobOptions.py
new file mode 100644
index 0000000000000000000000000000000000000000..c6e2898b4dd9e67738ec16ab773026fd012f745b
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/generatePRW_jobOptions.py
@@ -0,0 +1,23 @@
+
+theApp.EvtMax=-1 #says how many events to run over. Set to -1 for all events
+
+import AthenaPoolCnvSvc.ReadAthenaPool #sets up reading of POOL files (e.g. xAODs)
+from glob import glob
+svcMgr.EventSelector.InputCollections=glob( vars().get("FILES","/afs/cern.ch/atlas/project/PAT/xAODs/r5591/mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591/AOD.01494882._111853.pool.root.1").strip() )
+
+if len(svcMgr.EventSelector.InputCollections)==0:
+ print("WARNING >>>>>>>>>> generatePRW_jobOptions.py: NO INPUT FILES PROVIDED/FOUND FROM: %s ... this will produce a failure unless you are sending this job to the grid <<<<<<<<<<<<" % vars().get("FILES","/afs/cern.ch/atlas/project/PAT/xAODs/r5591/mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591/AOD.01494882._111853.pool.root.1").strip())
+
+
+algseq = CfgMgr.AthSequencer("AthAlgSeq") #gets the main AthSequencer
+algseq += CfgMgr.CP__PileupReweightingProvider(ConfigOutputStream="METADATA") #adds an instance of your alg to it
+
+include("AthAnalysisBaseComps/SuppressLogging.py") #Optional include to suppress as much athena output as possible
+
+svcMgr += CfgMgr.THistSvc()
+svcMgr.THistSvc.Output += ["METADATA DATAFILE='my.prw.root' OPT='RECREATE'"]
+
+#use on the grid like this:
+#pathena PileupReweighting/generatePRW_jobOptions.py --inDS="etc/,etc/,etc/" --outDS="user.whatever.myprw/" --extOutFile="auto.prw.root"
+
+
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/makeWeightTree.C b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/makeWeightTree.C
new file mode 100644
index 0000000000000000000000000000000000000000..d5afc46e43a1ad1e349e053c10fe020e13f90254
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/share/makeWeightTree.C
@@ -0,0 +1,37 @@
+
+///***************************************************************************************
+///PRW Weight Tree Making and Usage:
+///---------------------------------------------------------------------------------------
+///This is a small script you can run to generate a 'prw weight tree'
+///Modify the properties below to set up an instance of a prw tool as you would if you
+///processing all your MC in a usual event loop.
+///Then in the 'MakeWeightTree' line, specify all the channel numbers you want to generate
+///pileup weights for. Specify them as a comma delimited list in a string.
+///
+///Once you run this script (root makeWeightTree.C) you will get a prwTree.root file
+///To use it, you need to have made a branch in your main ntuple called 'PRWHash'
+///which contains the result of the GetPRWHash() method of the tool.
+///Then you just friend this tree to your main tree/chain, like this:
+///TChain c("MainTreeChain"); //add files to this
+///TChain* prwChain = new TChain("prwTree"); prwChain->Add("prwTree.root"); prwChain->BuildIndex("PRWHash");
+///c.AddFriend(prwChain);
+///Your main chain will then appear to have a 'PileupWeight' branch in it!!
+
+
+
+{
+ CP::PileupReweightingTool tool("tool");
+
+ std::vector<std::string> configFiles = {"/usera/will/testareas/ZdZd13TeV/ZdZdAnalysis/share/302073.prw.root","/usera/will/testareas/ZdZd13TeV/ZdZdAnalysis/share/361106_361108.r6630.prw.root"};
+ std::vector<std::string> lumiFiles = {"/usera/will/testareas/ZdZd13TeV/ZdZdAnalysis/share/data15_13TeV.periodA.None.lumicalc.gerl.root"};
+
+ tool.setProperty("ConfigFiles",configFiles);
+ tool.setProperty("LumicalcFiles",lumiFiles);
+ tool.setProperty("DefaultChannel",361106);
+
+
+ tool.initialize();
+ tool.expert()->MakeWeightTree("361106,361107,361108,123456","prwTree.root","PRWHash" /*change to name of branch in main tree*/ ,"PileupWeight" /*optional change to name of branch you want to be the weight*/);
+
+
+}
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/PileupReweightingProvider.h b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/PileupReweightingProvider.h
index 1196c13683f2a17294d6d04e28006396def18212..e13756df0caffff92d420ef68dca2b4818e88710 100644
--- a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/PileupReweightingProvider.h
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/PileupReweightingProvider.h
@@ -9,23 +9,83 @@
#include "AthenaBaseComps/AthAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
+#include "xAODEventInfo/EventInfoAuxContainer.h"
+
+#include "PATInterfaces/SystematicRegistry.h"
+
namespace CP {
class PileupReweightingProvider : public AthAlgorithm {
public:
- PileupReweightingProvider( const std::string& name, ISvcLocator* svcloc) : AthAlgorithm(name,svcloc),m_tool("PileupReweightingTool") {
+ PileupReweightingProvider( const std::string& name, ISvcLocator* svcloc) : AthAlgorithm(name,svcloc),m_tool("CP::PileupReweightingTool/auto") {
declareProperty("Tool",m_tool,"The configured PileupReweightingTool");
+ declareProperty("Input",m_inputKey="","Specify a specific EventInfo object");
+ declareProperty("Output",m_outputKey="","Specify an output EventInfo object. If differs from input, will create a clone of EventInfo and decorate that");
+ declareProperty("ConfigOutputStream",m_configStream="","Specify the stream to output config file to");
}
~PileupReweightingProvider() { }
- virtual StatusCode initialize() { CHECK( m_tool.retrieve() ); return StatusCode::SUCCESS; }
- virtual StatusCode execute() { CHECK( m_tool->execute() ); return StatusCode::SUCCESS; }
+ virtual StatusCode initialize() {
+ CHECK( m_tool.retrieve() );
+
+ IProperty* myTool = dynamic_cast<IProperty*>(&*m_tool);
+ CHECK( myTool->setProperty("ConfigOutputStream",m_configStream) );
+
+ //get the list of systematics
+ CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance();
+ m_allSysts = registry.recommendedSystematics();
+
+ if(m_configStream!="") ATH_MSG_INFO("Now running config file making .... please be patient! ... ");
+
+ return StatusCode::SUCCESS;
+ }
+ virtual StatusCode execute() {
+ const xAOD::EventInfo* evtInfo =0;
+
+ if(m_inputKey.length()>0) ATH_CHECK(evtStore()->retrieve(evtInfo,m_inputKey));
+ else {
+ #ifdef XAOD_STANDALONE
+ ATH_CHECK( evtStore()->retrieve(evtInfo,"") ); //apparently TEvent can't do keyless retrieves!!??
+ #else
+ ATH_CHECK(evtStore()->retrieve(evtInfo));
+ #endif
+ }
+
+ //do we need to make a copy??
+ if(m_inputKey!=m_outputKey && m_outputKey!="") {
+ xAOD::EventInfo* evtInfoCopy = new xAOD::EventInfo( *evtInfo );
+ xAOD::EventInfoAuxContainer* aux = new xAOD::EventInfoAuxContainer;
+ evtInfoCopy->setStore(aux);
+ ATH_CHECK( evtStore()->record( evtInfoCopy , m_outputKey ) );
+ ATH_CHECK( evtStore()->record( aux , m_outputKey+"Aux." ) );
+ evtInfo = evtInfoCopy;
+ }
+
+ CHECK( m_tool->apply(*evtInfo) );
+
+
+ //here's an example of systematic variations
+ for(auto& syst : m_allSysts) {
+ ATH_MSG_VERBOSE("Doing systematic : " << syst.name());
+ if(! m_tool->isAffectedBySystematic( syst )) continue;
+ CP::SystematicSet tmp; tmp.insert( syst );
+ if( m_tool->applySystematicVariation( tmp ) != CP::SystematicCode::Ok ) continue;
+ CHECK( m_tool->apply(*evtInfo) );
+ }
+ //make sure we leave the tool in the nominal mode
+ if( m_tool->applySystematicVariation( CP::SystematicSet() ) != CP::SystematicCode::Ok ) return StatusCode::FAILURE;
+
+ return StatusCode::SUCCESS;
+ }
private:
ToolHandle<IPileupReweightingTool> m_tool;
+ std::string m_inputKey,m_outputKey,m_configStream;
+
+ CP::SystematicSet m_allSysts;
};
diff --git a/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/testPRW.C b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/testPRW.C
new file mode 100644
index 0000000000000000000000000000000000000000..81f897fb2d63dc85b3a91b1f849bb725adf7acd6
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PileupReweighting/src/testPRW.C
@@ -0,0 +1,184 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "PileupReweighting/TPileupReweighting.h"
+#include "TTree.h"
+
+void testValue(double v1, double v2) {
+ if( (!v1 && v2) || (!v2 && v1) ) throw std::runtime_error("zero-non-zero mismatch");
+ if( fabs(v1-v2)/v1 > 0.001 ) throw std::runtime_error("mismatch > 0.1% in values");
+}
+
+
+int main() {
+ //create a 'dummy' lumicalc file, with only a few runs, with a few lumiblocks
+ TFile l("dummy.None.lumicalc.root","RECREATE");
+ TTree *t = new TTree("LumiMetaData","LumiMetaData");
+ UInt_t runNbr=0;Float_t ps1=1;Float_t ps2=1; Float_t ps3=1;UInt_t lbn=0;Float_t intLumi=0;Float_t mu=0;
+ t->Branch("RunNbr",&runNbr);t->Branch("L1Presc",&ps1);t->Branch("L2Presc",&ps2);t->Branch("L3Presc",&ps3);t->Branch("LBStart",&lbn);t->Branch("IntLumi",&intLumi);t->Branch("AvergeInteractionPerXing",&mu);
+
+ runNbr=1;lbn=1;intLumi=5;mu=0.5;t->Fill();
+ runNbr=1;lbn=2;intLumi=10;mu=1.5;t->Fill();
+ runNbr=2;lbn=1;intLumi=3;mu=1.5;t->Fill();
+ runNbr=3;lbn=1;intLumi=7;mu=1.5;t->Fill();
+ runNbr=3;lbn=2;intLumi=2;mu=2.5;t->Fill(); //will be unrepresented data in the setup in channel 2000
+
+ t->Write();delete t;
+ l.Close();
+
+ //make some trigger lumicalcs... only the prescale numbers are important!
+ TFile l2("dummy.TriggerA.lumicalc.root","RECREATE");
+ t = new TTree("LumiMetaData","LumiMetaData");
+ t->Branch("RunNbr",&runNbr);t->Branch("L1Presc",&ps1);t->Branch("L2Presc",&ps2);t->Branch("L3Presc",&ps3);t->Branch("LBStart",&lbn);t->Branch("IntLumi",&intLumi);t->Branch("AvergeInteractionPerXing",&mu);
+
+ runNbr=1;lbn=1;ps1=2;mu=0.5;t->Fill();
+ runNbr=1;lbn=2;ps1=3;mu=1.5;t->Fill();
+ runNbr=2;lbn=1;ps1=3;mu=1.5;t->Fill();
+ runNbr=3;lbn=1;ps1=4;mu=1.5;t->Fill();
+ runNbr=3;lbn=2;ps1=2;mu=2.5;t->Fill();
+
+ t->Write();delete t;
+ l2.Close();
+
+ TFile l3("dummy.TriggerB.lumicalc.root","RECREATE");
+ t = new TTree("LumiMetaData","LumiMetaData");
+ t->Branch("RunNbr",&runNbr);t->Branch("L1Presc",&ps1);t->Branch("L2Presc",&ps2);t->Branch("L3Presc",&ps3);t->Branch("LBStart",&lbn);t->Branch("IntLumi",&intLumi);t->Branch("AvergeInteractionPerXing",&mu);
+
+ runNbr=1;lbn=1;ps1=4;mu=0.5;t->Fill();
+ runNbr=1;lbn=2;ps1=6;mu=1.5;t->Fill();
+ runNbr=2;lbn=1;ps1=6;mu=1.5;t->Fill();
+ runNbr=3;lbn=1;ps1=8;mu=1.5;t->Fill();
+ runNbr=3;lbn=2;ps1=4;mu=2.5;t->Fill();
+
+ t->Write();delete t;
+ l3.Close();
+
+ //now create a config file ... this test config file generating
+
+ CP::TPileupReweighting g;
+ g.AddPeriod(100,1,2);
+ g.AddPeriod(101,3,3);
+ g.SetUniformBinning(10,0,10);
+ g.Initialize();
+
+ g.Fill(100,2002,1,0.5);
+ g.Fill(100,2002,1,1.5);
+ g.Fill(101,2002,1,1.5);
+ g.Fill(101,2002,1,2.5);
+ g.Fill(101,2002,1,3.5);
+
+ g.WriteToFile("dummy1.prw.root");
+
+ CP::TPileupReweighting* genConfig = new CP::TPileupReweighting("genConfig");
+ genConfig->AddPeriod(100,1,2);
+ genConfig->AddPeriod(101,3,3);
+ genConfig->SetUniformBinning(10,0,10);
+ genConfig->Initialize();
+
+ genConfig->Fill(100,2000/*channel number*/,1/*evt_weight*/,0.5);
+ genConfig->Fill(100,2000,1,1.5);
+ genConfig->Fill(101,2000,1,1.5);genConfig->Fill(101,2000,1,1.5);genConfig->Fill(101,2000,1,1.5); //three times stats in this bin
+
+ genConfig->Fill(100,2001,0.5,1.5);
+ genConfig->Fill(101,2001,2,1.5);
+ genConfig->Fill(101,2001,1,2.5);
+
+ genConfig->WriteToFile("dummy2.prw.root");
+
+ //using dummy1 alone should be ok, no unrepresented data...
+ CP::TPileupReweighting prw1;
+ prw1.AddConfigFile("dummy1.prw.root");
+ prw1.AddLumiCalcFile("dummy.None.lumicalc.root");
+ prw1.AddLumiCalcFile("dummy.TriggerA.lumicalc.root","TriggerA");
+ prw1.AddLumiCalcFile("dummy.TriggerB.lumicalc.root","TriggerB");
+ //prw1.EnableDebugging(1);
+ prw1.Initialize();
+
+ std::cout << "prw1 Integrated lumi = " << prw1.GetIntegratedLumi() << " (expected=2.7e-5) " << std::endl; testValue(prw1.GetIntegratedLumi(),2.7e-5);
+ std::cout << "prw1 periodWeights : " << prw1.GetPeriodWeight(100,2002) << " (expected=1.6666) " << prw1.GetPeriodWeight(101,2002) << " (expected=0.5555)" << std::endl; testValue(prw1.GetPeriodWeight(100,2002),1.6666);
+ std::cout << "prw1 primaryWeights : " << prw1.GetPrimaryWeight(100,2002,0.5) << " (expected=0.5555) " << prw1.GetPrimaryWeight(100,2002,1.5) << " (expected=1.4444) " << prw1.GetPrimaryWeight(101,2002,0.5) << " (expected=inf .. no MC) " << prw1.GetPrimaryWeight(101,2002,1.5) << " (expected=2.3333) " << prw1.GetPrimaryWeight(101,2002,2.5) << " (expected=0.6666) " << prw1.GetPrimaryWeight(101,2002,3.5) << " (expected=0.) " << std::endl;
+ testValue(prw1.GetPrimaryWeight(100,2002,0.5),0.5555);
+ testValue(prw1.GetPrimaryWeight(100,2002,1.5),1.4444);
+ testValue(prw1.GetPrimaryWeight(101,2002,1.5),2.3333);
+ testValue(prw1.GetPrimaryWeight(101,2002,2.5),0.6666);
+ testValue(prw1.GetPrimaryWeight(101,2002,3.5),0);
+
+ //test the random run number and lumiblock number generating
+ int num1(0); int lnum1(0);
+ for(int i=0;i<1000000;i++) {
+ if(prw1.GetRandomRunNumber(100)==1) { //returns 1 or 2 in this test case .. should give 1 about 15/18ths of the time
+ num1++;
+ if(prw1.GetRandomLumiBlockNumber(1)==1) lnum1++;
+ }
+ }
+ std::cout << "prw1 fraction of run1 = " << float(num1)/1000000 << " (expected=0.833333) " << " fraction of these with lbn1 = " << float(lnum1)/num1 << " (expected=0.333333) " << std::endl;
+ testValue(float(num1)/1000000,0.8333333);
+ testValue(float(lnum1)/num1,0.333333);
+
+ //test the DataWeights functionality .. first test the IntegratedLumi ...
+ std::cout << "prw1 TriggerA lumi = " << prw1.GetIntegratedLumi("TriggerA") << " (expected=9.58333e-6) " << std::endl; testValue(prw1.GetIntegratedLumi("TriggerA"),9.58333e-6);
+ std::cout << "prw1 TriggerB lumi = " << prw1.GetIntegratedLumi("TriggerB") << " (expected=4.79166e-6) " << std::endl; testValue(prw1.GetIntegratedLumi("TriggerB"),4.79166e-6);
+ std::cout << "prw1 TriggerA&&TriggerB lumi = " << prw1.GetIntegratedLumi("TriggerA&&TriggerB") << " (expected=1.81597e-6) " << std::endl;testValue(prw1.GetIntegratedLumi("TriggerA&&TriggerB"),1.81597e-6);
+
+ //now check the dataweight .. remember the runNumber is converted into a periodNumber (i.e. the weights aren't unique to a runNumber, they are unique to the periodNumber)
+ std::cout << "prw TriggerA DataWeights (mu independent): " << prw1.GetDataWeight(1,"TriggerA") << " (expected=2.63415) " << prw1.GetDataWeight(2,"TriggerA") << " (expected=2.63415) " << prw1.GetDataWeight(3,"TriggerA") << " (expected=3.2727)" << std::endl;
+ std::cout << "prw TriggerA DataWeights (mu dependent): " << prw1.GetDataWeight(1,"TriggerA",0.5) << " (expected=2) " << prw1.GetDataWeight(1,"TriggerA",1.5) << " (expected=3) " << prw1.GetDataWeight(2,"TriggerA",1.5) << " (expected=3.0000) " << prw1.GetDataWeight(3,"TriggerA",1.5) << " (expected=4) " << prw1.GetDataWeight(3,"TriggerA",2.5) << " (expected=2) " << prw1.GetDataWeight(3,"TriggerA",3.5) << " (expected=nan) " << std::endl;
+ testValue(prw1.GetDataWeight(1,"TriggerA"),2.63415);
+ testValue(prw1.GetDataWeight(2,"TriggerA"),2.63415);
+ testValue(prw1.GetDataWeight(3,"TriggerA"),3.2727);
+ testValue(prw1.GetDataWeight(1,"TriggerA",0.5),2.);
+ testValue(prw1.GetDataWeight(2,"TriggerA",1.5),3.);
+ testValue(prw1.GetDataWeight(3,"TriggerA",1.5),4);
+ testValue(prw1.GetDataWeight(3,"TriggerA",2.5),2);
+ //testValue(prw1.GetDataWeight(3,"TriggerA",3.5),nan);
+
+ try {
+ CP::TPileupReweighting t;
+ t.AddConfigFile("dummy2.prw.root");
+ t.AddLumiCalcFile("dummy.None.lumicalc.root");
+ t.Initialize();
+ } catch(...) {
+ std::cout << "exception thrown as expected" << std::endl;
+ }
+
+ CP::TPileupReweighting prw2;
+ prw2.AddConfigFile("dummy2.prw.root");
+ prw2.AddLumiCalcFile("dummy.None.lumicalc.root");
+ prw2.SetUnrepresentedDataAction(2,1.);
+ prw2.Initialize();
+
+ std::cout << "prw2 channel2000 combined weights: " << prw2.GetCombinedWeight(100,2000,0.5) << " (expected 0.9259) " << prw2.GetCombinedWeight(100,2000,1.5) << " (expected 2.40741) " << prw2.GetCombinedWeight(101,2000,1.5) << " (expected=0.432099) " << std::endl;
+ std::cout << "prw2 channel2000 'sum of weights': " << prw2.GetCombinedWeight(100,2000,0.5)+prw2.GetCombinedWeight(100,2000,1.5)+3*prw2.GetCombinedWeight(101,2000,1.5) << " from 5 events (sum of evtweights=5)" << std::endl;
+ std::cout << "prw2 channel2000 unrepresented weights: " << prw2.GetUnrepresentedDataWeight(100,2000) << " " << prw2.GetUnrepresentedDataWeight(101,2000) << std::endl;
+ std::cout << "prw2 channel2000 combined*unrepresented weight: " << prw2.GetCombinedWeight(100,2000,0.5)*prw2.GetUnrepresentedDataWeight(100,2000)+prw2.GetCombinedWeight(100,2000,1.5)*prw2.GetUnrepresentedDataWeight(100,2000)+3*prw2.GetCombinedWeight(101,2000,1.5)*prw2.GetUnrepresentedDataWeight(101,2000) << " from 5 events (sum of evtweights=5)" << std::endl;
+
+ std::cout << "prw2 channel2001 combined weights: " << prw2.GetCombinedWeight(100,2001,1.5) << " (expected 3.37037) " << prw2.GetCombinedWeight(101,2001,1.5) << " (expected 0.45370) " << prw2.GetCombinedWeight(101,2001,2.5) << " (expected=0.259259) " << std::endl;
+ std::cout << "prw2 channel2001 'sum of weights': " << 0.5*prw2.GetCombinedWeight(100,2001,1.5)+2.*prw2.GetCombinedWeight(101,2001,1.5)+1.*prw2.GetCombinedWeight(101,2001,2.5) << " from 3 events (of various weights, sum of evtweights=3.5)" << std::endl;
+ std::cout << "prw2 channel2001 unrepresented weight: " << prw2.GetUnrepresentedDataWeight(100,2001) << " " << prw2.GetUnrepresentedDataWeight(101,2001) << std::endl;
+ std::cout << "prw2 channel2001 combined*unrepresented weight: " << 0.5*prw2.GetCombinedWeight(100,2001,1.5)*prw2.GetUnrepresentedDataWeight(100,2001)+2.*prw2.GetCombinedWeight(101,2001,1.5)*prw2.GetUnrepresentedDataWeight(101,2001)+1.*prw2.GetCombinedWeight(101,2001,2.5)*prw2.GetUnrepresentedDataWeight(101,2001) << " from 3 events (of various weights, sum of evtweights=3.5)" << std::endl;
+
+ //repeat with action=3 ...
+ CP::TPileupReweighting prw3;
+ prw3.AddConfigFile("dummy2.prw.root");
+ prw3.AddLumiCalcFile("dummy.None.lumicalc.root");
+ prw3.SetUnrepresentedDataAction(3,1.);
+ prw3.Initialize();
+
+ std::cout << "prw3 channel2000 combined weights: " << prw3.GetCombinedWeight(100,2000,0.5) << " (expected 0.0) " << prw3.GetCombinedWeight(100,2000,1.5) << " (expected 3.33333) " << prw3.GetCombinedWeight(101,2000,1.5) << " (expected=0.555555) " << std::endl;
+ std::cout << "prw3 channel2000 'sum of weights': " << prw3.GetCombinedWeight(100,2000,0.5)+prw3.GetCombinedWeight(100,2000,1.5)+3*prw3.GetCombinedWeight(101,2000,1.5) << " from 5 events (sum of evtweights=5)" << std::endl;
+
+ std::cout << "prw3 channel2001 combined weights: " << prw3.GetCombinedWeight(100,2001,1.5) << " (expected 4.66666) " << prw3.GetCombinedWeight(101,2001,1.5) << " (expected 0.58333) " << prw3.GetCombinedWeight(101,2001,2.5) << " (expected=0.259259) " << std::endl;
+ std::cout << "prw3 channel2001 'sum of weights': " << 0.5*prw3.GetCombinedWeight(100,2001,1.5)+2.*prw3.GetCombinedWeight(101,2001,1.5)+1.*prw3.GetCombinedWeight(101,2001,2.5) << " from 3 events (of various weights, sum of evtweights=3.5)" << std::endl;
+
+
+ //CP::TPileupReweighting t;
+
+ //t.AddLumiCalcFile("/usera/will/testareas/ZZ8TeV/ZZAnalysis/share/GRL_v61.None.lumicalc.root");
+
+
+
+
+ return 0;
+}
\ No newline at end of file