Utils/interface/PhotonFix.h

#ifndef PhotonFix_Defined_hh
#define PhotonFix_Defined_hh

//-------------------------------------------------------//
// Project:  PhotonFix
// Author:   Paul Dauncey (p.dauncey@imperial.ac.uk)
// Modified: 11/07/2011
// Admins:   Paul Dauncey (p.dauncey@imperial.ac.uk)
//           Matt Kenzie (matthew.william.kenzie@cern.ch)
//-------------------------------------------------------//

/*
  Does post-reco fixes to ECAL photon energy and estimates resolution.
  This can run outside of the usual CMS software framework but requires 
  access to a file 'EcalGaps.dat' which must be in the same directory as 
  that used to run.

  To run within CMSSW use PhotonFixCMS.h (which can access the geometry 
  directly - go to "RecoEcal/EgammaCoreTools/plugins/PhotonFixCMS.h"
  for details.

  Before instantiating any objects of PhotonFix, the constants must be
  initialised in the first event using
    PhotonFix::initialise("3_8");
  
  The string gives the reco version used. Valid strings are 
  "3_8", "3_11", "4_2" and "Nominal", where the latter gives no correction 
  to the energy and a nominal resolution value. There is also "4_2e" which 
  provides corrections for electrons which are reconstructed as photons (to
  aid with testing the performance of these corrections in data).

  Make objects using
    PhotonFix a(energy,eta,phi,r9);
  where energy is the photon energy, eta and phi are the ECAL
  cluster positions (NB from the Supercluster object, _not_ the
  Photon object, as the latter gives eta and phi directions,
  not positions), and r9 is the R9 value of the SC.

  Get the corrected energy using
    a.fixedEnergy();
  and the resolution using
    a.sigmaEnergy();

*/

#include <iostream>
#include <string>

class PhotonFix {
 public:
  //PhotonFix(double e, double eta, double phi, double r9);

  PhotonFix() : _onePi(acos(-1.0)), _twoPi(2.0*acos(-1.0)), _initialised(false) {}
  // Must be called before instantiating any PhotonFix objects
  bool initialise(const std::string &s="Nominal", const std::string &infile="PhotonFix.dat");
  bool initialised() ;
  
  void setup(double e, double eta, double phi, double r9);
  
  bool isbarrel() const { return !_be; }
  
  // Corrected energy and sigma
  double fixedEnergy() const;
  double sigmaEnergy() const;
  
  // Input values
  double rawEnergy() const;
  double eta() const;
  double phi() const;
  double r9() const;

  // Derived EB crystal, submodule and module relative coordinates
  double etaC() const;
  double etaS() const;
  double etaM() const;

  double phiC() const;
  double phiS() const;
  double phiM() const;

  // Derived EE zeta, crystal, subcrystal and D-module relative coordinates
  double xZ() const;
  double xC() const;
  double xS() const;
  double xM() const;

  double yZ() const;
  double yC() const;
  double yS() const;
  double yM() const;

  // Return arrays containing positions of ecal gaps
  void barrelCGap(unsigned i, unsigned j, unsigned k, double c);
  void barrelSGap(unsigned i, unsigned j, unsigned k, double c);
  void barrelMGap(unsigned i, unsigned j, unsigned k, double c);
  void endcapCrystal(unsigned i, unsigned j, bool c); 
  void endcapCGap(unsigned i, unsigned j, unsigned k, double c);
  void endcapSGap(unsigned i, unsigned j, unsigned k, double c);
  void endcapMGap(unsigned i, unsigned j, unsigned k, double c);
  
  void print() const;

  // Input and output the fit parameters
  void setParameters(unsigned be, unsigned hl, const double *p);
  void getParameters(unsigned be, unsigned hl, double *p);

  void dumpParameters(std::ostream &o);
  void printParameters(std::ostream &o);

  // Utility functions
  double GetaPhi(double f0, double f1) const;
  double asinh(double s) const;  
  void dumpGaps(std::ostream &o);
  
 private:


  // Used by above; do not call directly
  bool initialiseParameters(const std::string &s);
  bool initialiseGeometry(const std::string &s, const std::string &infile);

  
  // Utility functions
  double dPhi(double f0, double f1) const;
  double aPhi(double f0, double f1) const;

  double expCorrection(double a, const double *p) const;
  double gausCorrection(double a, const double *p) const;

  // Actual data for each instantiated object
  unsigned _be,_hl;
  double _e,_eta,_phi,_r9;
  double _aC,_aS,_aM,_bC,_bS,_bM;
  
  // Constants
  const double _onePi;
  const double _twoPi;
  
  // Initialisation flag
  bool _initialised;
  
  // Parameters for fixes
  double _meanScale[2][2][4];
  double _meanAT[2][2][4];
  double _meanAC[2][2][4];
  double _meanAS[2][2][4];
  double _meanAM[2][2][4];
  double _meanBT[2][2][4];
  double _meanBC[2][2][4];
  double _meanBS[2][2][4];
  double _meanBM[2][2][4];
  double _meanR9[2][2][4];
  
  // Parameters for resolution
  double _sigmaScale[2][2][4];
  double _sigmaAT[2][2][4];
  double _sigmaAC[2][2][4];
  double _sigmaAS[2][2][4];
  double _sigmaAM[2][2][4];
  double _sigmaBT[2][2][4];
  double _sigmaBC[2][2][4];
  double _sigmaBS[2][2][4];
  double _sigmaBM[2][2][4];
  double _sigmaR9[2][2][4];
  
  // EB gap positions
  double _barrelCGap[169][360][2];
  double _barrelSGap[33][180][2];
  double _barrelMGap[7][18][2];
  
  // EE crystal existence and gap positions
  bool   _endcapCrystal[100][100];
  double _endcapCGap[2][7080][2];
  double _endcapSGap[2][264][2];
  double _endcapMGap[2][1][2];

};

#endif
Revision:	1.2
Committed:	Thu Sep 8 15:51:23 2011 UTC (13 years, 7 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_029c, Mit_029b, Mit_029a, Mit_028a, Mit_028, Mit_027, Mit_027a, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_025pre1, HEAD
Changes since 1.1:	+68 -57 lines
Log Message:	Add tool and required changes for photon energy regression
#	User	Rev	Content
1	bendavid	1.1	#ifndef PhotonFix_Defined_hh
2			#define PhotonFix_Defined_hh
3
4			//-------------------------------------------------------//
5			// Project: PhotonFix
6			// Author: Paul Dauncey (p.dauncey@imperial.ac.uk)
7			// Modified: 11/07/2011
8			// Admins: Paul Dauncey (p.dauncey@imperial.ac.uk)
9			// Matt Kenzie (matthew.william.kenzie@cern.ch)
10			//-------------------------------------------------------//
11
12			/*
13			Does post-reco fixes to ECAL photon energy and estimates resolution.
14			This can run outside of the usual CMS software framework but requires
15	bendavid	1.2	access to a file 'EcalGaps.dat' which must be in the same directory as
16	bendavid	1.1	that used to run.
17
18			To run within CMSSW use PhotonFixCMS.h (which can access the geometry
19			directly - go to "RecoEcal/EgammaCoreTools/plugins/PhotonFixCMS.h"
20			for details.
21
22			Before instantiating any objects of PhotonFix, the constants must be
23			initialised in the first event using
24			PhotonFix::initialise("3_8");
25
26			The string gives the reco version used. Valid strings are
27			"3_8", "3_11", "4_2" and "Nominal", where the latter gives no correction
28	bendavid	1.2	to the energy and a nominal resolution value. There is also "4_2e" which
29			provides corrections for electrons which are reconstructed as photons (to
30			aid with testing the performance of these corrections in data).
31	bendavid	1.1
32			Make objects using
33			PhotonFix a(energy,eta,phi,r9);
34			where energy is the photon energy, eta and phi are the ECAL
35			cluster positions (NB from the Supercluster object, _not_ the
36			Photon object, as the latter gives eta and phi directions,
37			not positions), and r9 is the R9 value of the SC.
38
39			Get the corrected energy using
40			a.fixedEnergy();
41			and the resolution using
42			a.sigmaEnergy();
43
44			*/
45
46			#include <iostream>
47			#include <string>
48
49			class PhotonFix {
50			public:
51	bendavid	1.2	//PhotonFix(double e, double eta, double phi, double r9);
52	bendavid	1.1
53	bendavid	1.2	PhotonFix() : _onePi(acos(-1.0)), _twoPi(2.0*acos(-1.0)), _initialised(false) {}
54	bendavid	1.1	// Must be called before instantiating any PhotonFix objects
55	bendavid	1.2	bool initialise(const std::string &s="Nominal", const std::string &infile="PhotonFix.dat");
56			bool initialised() ;
57
58			void setup(double e, double eta, double phi, double r9);
59
60			bool isbarrel() const { return !_be; }
61	bendavid	1.1
62			// Corrected energy and sigma
63			double fixedEnergy() const;
64			double sigmaEnergy() const;
65
66			// Input values
67			double rawEnergy() const;
68			double eta() const;
69			double phi() const;
70			double r9() const;
71
72			// Derived EB crystal, submodule and module relative coordinates
73			double etaC() const;
74			double etaS() const;
75			double etaM() const;
76
77			double phiC() const;
78			double phiS() const;
79			double phiM() const;
80
81			// Derived EE zeta, crystal, subcrystal and D-module relative coordinates
82			double xZ() const;
83			double xC() const;
84			double xS() const;
85			double xM() const;
86
87			double yZ() const;
88			double yC() const;
89			double yS() const;
90			double yM() const;
91
92			// Return arrays containing positions of ecal gaps
93	bendavid	1.2	void barrelCGap(unsigned i, unsigned j, unsigned k, double c);
94			void barrelSGap(unsigned i, unsigned j, unsigned k, double c);
95			void barrelMGap(unsigned i, unsigned j, unsigned k, double c);
96			void endcapCrystal(unsigned i, unsigned j, bool c);
97			void endcapCGap(unsigned i, unsigned j, unsigned k, double c);
98			void endcapSGap(unsigned i, unsigned j, unsigned k, double c);
99			void endcapMGap(unsigned i, unsigned j, unsigned k, double c);
100	bendavid	1.1
101			void print() const;
102
103			// Input and output the fit parameters
104	bendavid	1.2	void setParameters(unsigned be, unsigned hl, const double *p);
105			void getParameters(unsigned be, unsigned hl, double *p);
106	bendavid	1.1
107	bendavid	1.2	void dumpParameters(std::ostream &o);
108			void printParameters(std::ostream &o);
109	bendavid	1.1
110			// Utility functions
111	bendavid	1.2	double GetaPhi(double f0, double f1) const;
112			double asinh(double s) const;
113			void dumpGaps(std::ostream &o);
114	bendavid	1.1
115			private:
116
117	bendavid	1.2
118			// Used by above; do not call directly
119			bool initialiseParameters(const std::string &s);
120			bool initialiseGeometry(const std::string &s, const std::string &infile);
121
122
123
124
125	bendavid	1.1	// Utility functions
126	bendavid	1.2	double dPhi(double f0, double f1) const;
127			double aPhi(double f0, double f1) const;
128	bendavid	1.1
129	bendavid	1.2	double expCorrection(double a, const double *p) const;
130			double gausCorrection(double a, const double *p) const;
131	bendavid	1.1
132			// Actual data for each instantiated object
133			unsigned _be,_hl;
134			double _e,_eta,_phi,_r9;
135			double _aC,_aS,_aM,_bC,_bS,_bM;
136
137			// Constants
138	bendavid	1.2	const double _onePi;
139			const double _twoPi;
140	bendavid	1.1
141			// Initialisation flag
142	bendavid	1.2	bool _initialised;
143	bendavid	1.1
144			// Parameters for fixes
145	bendavid	1.2	double _meanScale[2][2][4];
146			double _meanAT[2][2][4];
147			double _meanAC[2][2][4];
148			double _meanAS[2][2][4];
149			double _meanAM[2][2][4];
150			double _meanBT[2][2][4];
151			double _meanBC[2][2][4];
152			double _meanBS[2][2][4];
153			double _meanBM[2][2][4];
154			double _meanR9[2][2][4];
155	bendavid	1.1
156			// Parameters for resolution
157	bendavid	1.2	double _sigmaScale[2][2][4];
158			double _sigmaAT[2][2][4];
159			double _sigmaAC[2][2][4];
160			double _sigmaAS[2][2][4];
161			double _sigmaAM[2][2][4];
162			double _sigmaBT[2][2][4];
163			double _sigmaBC[2][2][4];
164			double _sigmaBS[2][2][4];
165			double _sigmaBM[2][2][4];
166			double _sigmaR9[2][2][4];
167	bendavid	1.1
168			// EB gap positions
169	bendavid	1.2	double _barrelCGap[169][360][2];
170			double _barrelSGap[33][180][2];
171			double _barrelMGap[7][18][2];
172	bendavid	1.1
173			// EE crystal existence and gap positions
174	bendavid	1.2	bool _endcapCrystal[100][100];
175			double _endcapCGap[2][7080][2];
176			double _endcapSGap[2][264][2];
177			double _endcapMGap[2][1][2];
178	bendavid	1.1
179			};
180
181			#endif