PatAlgos/plugins/ObjectEnergyScale.h

//
// $Id: ObjectEnergyScale.h,v 1.1 2008/03/06 09:23:10 llista Exp $
//

#ifndef PhysicsTools_PatAlgos_ObjectEnergyScale_h
#define PhysicsTools_PatAlgos_ObjectEnergyScale_h

/**
  \class    pat::ObjectEnergyScale ObjectEnergyScale.h "PhysicsTools/PatAlgos/interface/ObjectEnergyScale.h"
  \brief    Energy scale shifting and smearing module

   This class provides energy scale shifting & smearing to objects with
   resolutions for systematic error studies.
   A detailed documentation is found in
     PhysicsTools/PatAlgos/data/ObjectEnergyScale.cfi

  \author   Volker Adler
  \version  $Id: ObjectEnergyScale.h,v 1.1 2008/03/06 09:23:10 llista Exp $
*/


#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/InputTag.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "CLHEP/Random/RandGaussQ.h"


namespace pat {


  template<class T>
  class ObjectEnergyScale : public edm::EDProducer {

    public:

      explicit ObjectEnergyScale(const edm::ParameterSet& iConfig);
      ~ObjectEnergyScale();

    private:

      virtual void produce( edm::Event& iEvent, const edm::EventSetup& iSetup);

      float  getSmearing(T& object);
      void   setScale(T& object);

      edm::InputTag objects_;
      float         factor_,
                    shiftFactor_,
                    iniRes_,
                    worsenRes_;
      bool          useFixedMass_,
                    useDefaultIniRes_,
                    useIniResByFraction_,
                    useWorsenResByFactor_;

      CLHEP::RandGaussQ* gaussian_;

  };


}


template<class T>
pat::ObjectEnergyScale<T>::ObjectEnergyScale(const edm::ParameterSet& iConfig)
{
  objects_              = iConfig.getParameter<edm::InputTag>("scaledObject");
  useFixedMass_         = iConfig.getParameter<bool>         ("fixMass");
  shiftFactor_          = iConfig.getParameter<double>       ("shiftFactor");
  useDefaultIniRes_     = iConfig.getParameter<bool>         ("useDefaultInitialResolution");
  iniRes_               = iConfig.getParameter<double>       ("initialResolution");
  useIniResByFraction_  = iConfig.getParameter<bool>         ("initialResolutionByFraction");
  worsenRes_            = iConfig.getParameter<double>       ("worsenResolution");
  useWorsenResByFactor_ = iConfig.getParameter<bool>         ("worsenResolutionByFactor");

  edm::Service<edm::RandomNumberGenerator> rng;
  CLHEP::HepRandomEngine& engine = rng->getEngine();
  gaussian_ = new CLHEP::RandGaussQ(engine);

  produces<std::vector<T> >();
}


template<class T>
pat::ObjectEnergyScale<T>::~ObjectEnergyScale()
{
  delete gaussian_;
}


template<class T>
void pat::ObjectEnergyScale<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  edm::Handle<std::vector<T> > objectsHandle;
  iEvent.getByLabel(objects_, objectsHandle);
  std::vector<T> objects = *objectsHandle;
  std::auto_ptr<std::vector<T> > objectsVector(new std::vector<T>);
  objectsVector->reserve(objectsHandle->size());

  for ( unsigned int i = 0; i < objects.size(); i++ ) {
    factor_ = shiftFactor_ * ( objects[i].energy() > 0. ?
                               getSmearing(objects[i])  :
                               0.);
    setScale(objects[i]);
    objectsVector->push_back(objects[i]);
  }
  iEvent.put(objectsVector);
}


/// Returns a smearing factor which is multiplied to the initial value then to get it smeared,
/// sets initial resolution to resolution provided by input object if required
/// and converts the 'worsenResolution' parameter to protect from meaningless final resolution values.
template<class T>
float pat::ObjectEnergyScale<T>::getSmearing(T& object)
{
  // overwrite config file parameter 'initialResolution' if required
  if ( useDefaultIniRes_ ) {
    // get initial resolution from input object (and calculate relative initial resolution from absolute value)
    iniRes_ = (1. / sin(object.theta()) * object.resolutionEt() - object.et() * cos(object.theta()) / pow(sin(object.theta()),2) * object.resolutionTheta()) / object.energy(); // conversion of resEt and resTheta into energy resolution
  } else if ( ! useIniResByFraction_ ) {
    // calculate relative initial resolution from absolute value
    iniRes_ = iniRes_ / object.energy();
  }
  // Is 'worsenResolution' a factor or a summand?
  float finalRes = useWorsenResByFactor_                            ?
                    (1.+fabs(1.-fabs(worsenRes_)))   * fabs(iniRes_) :
                    fabs(worsenRes_)/object.energy() + fabs(iniRes_); // conversion as protection from "finalRes_<iniRes_"
  // return smearing factor
  return std::max( gaussian_->fire(1., sqrt(pow(finalRes,2)-pow(iniRes_,2))), 0. ); // protection from negative smearing factors
}


/// Mutliplies the final factor (consisting of shifting and smearing factors) to the object's 4-vector
/// and takes care of preserved masses.
template<class T>
void pat::ObjectEnergyScale<T>::setScale(T& object)
{
  if ( factor_ < 0. ) {
    factor_ = 0.;
  }
  // calculate the momentum factor for fixed or not fixed mass
  float factorMomentum = useFixedMass_ && object.p() > 0.                                   ?
                          sqrt(pow(factor_*object.energy(),2)-object.massSqr()) / object.p() :
                          factor_;
  // set shifted & smeared new 4-vector
  object.setP4(reco::Particle::LorentzVector(factorMomentum*object.px(),
                                             factorMomentum*object.py(),
                                             factorMomentum*object.pz(),
                                             factor_       *object.energy()));
}


#endif
Revision:	1.1
Committed:	Tue Oct 20 17:15:14 2009 UTC (15 years, 6 months ago) by dgele
Content type:	text/plain
Branch:	MAIN
Branch point for:	ANA
Log Message:	Initial revision
#	User	Rev	Content
1	dgele	1.1	//
2			// $Id: ObjectEnergyScale.h,v 1.1 2008/03/06 09:23:10 llista Exp $
3			//
4
5			#ifndef PhysicsTools_PatAlgos_ObjectEnergyScale_h
6			#define PhysicsTools_PatAlgos_ObjectEnergyScale_h
7
8			/**
9			\class pat::ObjectEnergyScale ObjectEnergyScale.h "PhysicsTools/PatAlgos/interface/ObjectEnergyScale.h"
10			\brief Energy scale shifting and smearing module
11
12			This class provides energy scale shifting & smearing to objects with
13			resolutions for systematic error studies.
14			A detailed documentation is found in
15			PhysicsTools/PatAlgos/data/ObjectEnergyScale.cfi
16
17			\author Volker Adler
18			\version $Id: ObjectEnergyScale.h,v 1.1 2008/03/06 09:23:10 llista Exp $
19			*/
20
21
22			#include "FWCore/Framework/interface/Frameworkfwd.h"
23			#include "FWCore/Framework/interface/EDProducer.h"
24			#include "FWCore/Framework/interface/Event.h"
25			#include "FWCore/ParameterSet/interface/ParameterSet.h"
26			#include "FWCore/ParameterSet/interface/InputTag.h"
27
28			#include "FWCore/ServiceRegistry/interface/Service.h"
29			#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
30			#include "CLHEP/Random/RandGaussQ.h"
31
32
33			namespace pat {
34
35
36			template<class T>
37			class ObjectEnergyScale : public edm::EDProducer {
38
39			public:
40
41			explicit ObjectEnergyScale(const edm::ParameterSet& iConfig);
42			~ObjectEnergyScale();
43
44			private:
45
46			virtual void produce( edm::Event& iEvent, const edm::EventSetup& iSetup);
47
48			float getSmearing(T& object);
49			void setScale(T& object);
50
51			edm::InputTag objects_;
52			float factor_,
53			shiftFactor_,
54			iniRes_,
55			worsenRes_;
56			bool useFixedMass_,
57			useDefaultIniRes_,
58			useIniResByFraction_,
59			useWorsenResByFactor_;
60
61			CLHEP::RandGaussQ* gaussian_;
62
63			};
64
65
66			}
67
68
69			template<class T>
70			pat::ObjectEnergyScale<T>::ObjectEnergyScale(const edm::ParameterSet& iConfig)
71			{
72			objects_ = iConfig.getParameter<edm::InputTag>("scaledObject");
73			useFixedMass_ = iConfig.getParameter<bool> ("fixMass");
74			shiftFactor_ = iConfig.getParameter<double> ("shiftFactor");
75			useDefaultIniRes_ = iConfig.getParameter<bool> ("useDefaultInitialResolution");
76			iniRes_ = iConfig.getParameter<double> ("initialResolution");
77			useIniResByFraction_ = iConfig.getParameter<bool> ("initialResolutionByFraction");
78			worsenRes_ = iConfig.getParameter<double> ("worsenResolution");
79			useWorsenResByFactor_ = iConfig.getParameter<bool> ("worsenResolutionByFactor");
80
81			edm::Service<edm::RandomNumberGenerator> rng;
82			CLHEP::HepRandomEngine& engine = rng->getEngine();
83			gaussian_ = new CLHEP::RandGaussQ(engine);
84
85			produces<std::vector<T> >();
86			}
87
88
89			template<class T>
90			pat::ObjectEnergyScale<T>::~ObjectEnergyScale()
91			{
92			delete gaussian_;
93			}
94
95
96			template<class T>
97			void pat::ObjectEnergyScale<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
98			{
99			edm::Handle<std::vector<T> > objectsHandle;
100			iEvent.getByLabel(objects_, objectsHandle);
101			std::vector<T> objects = *objectsHandle;
102			std::auto_ptr<std::vector<T> > objectsVector(new std::vector<T>);
103			objectsVector->reserve(objectsHandle->size());
104
105			for ( unsigned int i = 0; i < objects.size(); i++ ) {
106			factor_ = shiftFactor_ * ( objects[i].energy() > 0. ?
107			getSmearing(objects[i]) :
108			0.);
109			setScale(objects[i]);
110			objectsVector->push_back(objects[i]);
111			}
112			iEvent.put(objectsVector);
113			}
114
115
116			/// Returns a smearing factor which is multiplied to the initial value then to get it smeared,
117			/// sets initial resolution to resolution provided by input object if required
118			/// and converts the 'worsenResolution' parameter to protect from meaningless final resolution values.
119			template<class T>
120			float pat::ObjectEnergyScale<T>::getSmearing(T& object)
121			{
122			// overwrite config file parameter 'initialResolution' if required
123			if ( useDefaultIniRes_ ) {
124			// get initial resolution from input object (and calculate relative initial resolution from absolute value)
125			iniRes_ = (1. / sin(object.theta()) * object.resolutionEt() - object.et() * cos(object.theta()) / pow(sin(object.theta()),2) * object.resolutionTheta()) / object.energy(); // conversion of resEt and resTheta into energy resolution
126			} else if ( ! useIniResByFraction_ ) {
127			// calculate relative initial resolution from absolute value
128			iniRes_ = iniRes_ / object.energy();
129			}
130			// Is 'worsenResolution' a factor or a summand?
131			float finalRes = useWorsenResByFactor_ ?
132			(1.+fabs(1.-fabs(worsenRes_))) * fabs(iniRes_) :
133			fabs(worsenRes_)/object.energy() + fabs(iniRes_); // conversion as protection from "finalRes_<iniRes_"
134			// return smearing factor
135			return std::max( gaussian_->fire(1., sqrt(pow(finalRes,2)-pow(iniRes_,2))), 0. ); // protection from negative smearing factors
136			}
137
138
139			/// Mutliplies the final factor (consisting of shifting and smearing factors) to the object's 4-vector
140			/// and takes care of preserved masses.
141			template<class T>
142			void pat::ObjectEnergyScale<T>::setScale(T& object)
143			{
144			if ( factor_ < 0. ) {
145			factor_ = 0.;
146			}
147			// calculate the momentum factor for fixed or not fixed mass
148			float factorMomentum = useFixedMass_ && object.p() > 0. ?
149			sqrt(pow(factor_*object.energy(),2)-object.massSqr()) / object.p() :
150			factor_;
151			// set shifted & smeared new 4-vector
152			object.setP4(reco::Particle::LorentzVector(factorMomentum*object.px(),
153			factorMomentum*object.py(),
154			factorMomentum*object.pz(),
155			factor_ *object.energy()));
156			}
157
158
159			#endif