VHbbAnalysisCode/interface/LeptonInfoPlotSet.h

#ifndef trkupgradeanalysis_LeptonInfoPlotSet_h
#define trkupgradeanalysis_LeptonInfoPlotSet_h

#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"

// Forward declarations
class TH1F;
class TDirectory;
class TTree;
class VHbbEventAuxInfo;

namespace trkupgradeanalysis
{
        /** @brief An abstract class that incorporates some of the code used for both electron and muon plotsets.
         *
         * Note that this isn't intended to be used on its own, rather as a member for ElectronInfoPlotSet and
         * MuonInfoPlotSet classes to encapsulate common code between them.
         *
         * @author Mark Grimes (mark.grimes@bristol.ac.uk)
         * @date 02/Jul/2012
         */
        template<class T_Lepton>
        class LeptonInfoPlotSet
        {
        public:
                /** @brief Only constructor
                 *
                 * createIsolationTree - set whether to record the isolation values in a TTree or not. Default
                 * is off to save space.
                 */
                LeptonInfoPlotSet( bool createIsolationTree=false );
                void book( TDirectory* pDirectory );
                void fill( const T_Lepton& lepton, const VHbbEventAuxInfo* pAuxInfo=NULL );

                static float combinedRelativeIsolation( const T_Lepton& lepton );
                static float deltaBetaCorrectedIsolation( const T_Lepton& lepton, float deltaBetaFactor=-0.5 );
                static float rho25CorrectedIsolation( const T_Lepton& lepton, float rho25 );
                static float thirdRho25CorrectedIsolation( const T_Lepton& lepton, float rho25 );
                static float deltaBetaCorrectedIsolationNoZeroing( const T_Lepton& lepton, float deltaBetaFactor=-0.5 );
                static float rho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 );
                static float thirdRho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 );
        private:
                bool histogramHaveBeenBooked_;

                bool createIsolationTree_; ///< Whether or not to record isolation values in a TTree.
                TTree* pIsolationTree_; ///< This TTree is only created if createIsolationTree is set to true in the constructor.
                int numberOfPrimaryVertices_branch_; ///< All these are variables to hold the branch data
                float chargedIsolation_branch_;
                float photonIsolation_branch_;
                float neutralIsolation_branch_;
                float pileupIsolation_branch_;
                float rho25_branch_;
                float pT_branch_;

                TH1F* pEta_;
                TH1F* pPhi_;
                TH1F* pPt_;
                TH1F* pChargedIsolation_;
                TH1F* pPhotonIsolation_;
                TH1F* pNeutralIsolation_;
                TH1F* pPileupIsolation_;
                TH1F* pRelativeChargedIsolation_;
                TH1F* pRelativePhotonIsolation_;
                TH1F* pRelativeNeutralIsolation_;
                TH1F* pRelativePileupIsolation_;
                TH1F* pRelativeIsolation_;
                TH1F* pDeltaBetaCorrectedIsolation_;
                TH1F* pRho25CorrectedIsolation_;
                TH1F* pThirdRho25CorrectedIsolation_;
                TH1F* pDeltaBetaCorrectedIsolationNoZeroing_;
                TH1F* pRho25CorrectedIsolationNoZeroing_;
                TH1F* pThirdRho25CorrectedIsolationNoZeroing_;
        };

} // end of namespace trkupgradeanalysis


//----------------------------------------------------------
//
// Need the implementation here since it's templated
//
//----------------------------------------------------------
#include <stdexcept>
#include <TH1F.h>
#include <TTree.h>
#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"


template<class T_Lepton>
trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::LeptonInfoPlotSet( bool createIsolationTree )
        : histogramHaveBeenBooked_(false), createIsolationTree_(createIsolationTree), pIsolationTree_(NULL)
{
        // No operation besides the initialiser list.
}

template<class T_Lepton>
void trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::book( TDirectory* pDirectory )
{
        if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::LeptonInfoPlotSet::book() - histograms have already been booked" );

        //
        // Note that the root file which TDirectory is part of takes ownership of all
        // of these objects, so I don't need to (and shouldn't) delete them when I'm
        // finished.
        //


        pEta_=new TH1F( "eta", "Eta", 60, -3, 3 );
        pEta_->SetDirectory(pDirectory);

        pPhi_=new TH1F( "phi", "phi", 60, -3.15, 3.15 );
        pPhi_->SetDirectory(pDirectory);

        pPt_=new TH1F( "pT", "pT", 60, 0, 250 );
        pPt_->SetDirectory(pDirectory);

        pChargedIsolation_=new TH1F( "chargedIsolation", "chargedIsolation", 60, 0, 30 );
        pChargedIsolation_->SetDirectory(pDirectory);

        pPhotonIsolation_=new TH1F( "photonIsolation", "photonIsolation", 60, 0, 30 );
        pPhotonIsolation_->SetDirectory(pDirectory);

        pNeutralIsolation_=new TH1F( "neutralIsolation", "neutralIsolation", 60, 0, 30 );
        pNeutralIsolation_->SetDirectory(pDirectory);

        pPileupIsolation_=new TH1F( "pileupIsolation", "pileupIsolation", 60, 0, 30 );
        pPileupIsolation_->SetDirectory(pDirectory);

        pRelativeChargedIsolation_=new TH1F( "relativeChargedIsolation", "relativeChargedIsolation", 60, 0, 1 );
        pRelativeChargedIsolation_->SetDirectory(pDirectory);

        pRelativePhotonIsolation_=new TH1F( "relativePhotonIsolation", "relativePhotonIsolation", 60, 0, 1 );
        pRelativePhotonIsolation_->SetDirectory(pDirectory);

        pRelativeNeutralIsolation_=new TH1F( "relativeNeutralIsolation", "relativeNeutralIsolation", 60, 0, 1 );
        pRelativeNeutralIsolation_->SetDirectory(pDirectory);

        pRelativePileupIsolation_=new TH1F( "relativePileupIsolation", "relativePileupIsolation", 60, 0, 1 );
        pRelativePileupIsolation_->SetDirectory(pDirectory);

        pRelativeIsolation_=new TH1F( "relativeIsolation", "relativeIsolation", 120, 0, 2 );
        pRelativeIsolation_->SetDirectory(pDirectory);

        pDeltaBetaCorrectedIsolation_=new TH1F( "deltaBetaCorrectedIsolation", "deltaBetaCorrectedIsolation", 120, 0, 2 );
        pDeltaBetaCorrectedIsolation_->SetDirectory(pDirectory);

        pRho25CorrectedIsolation_=new TH1F( "rho25CorrectedIsolation", "rho25CorrectedIsolation", 120, 0, 2 );
        pRho25CorrectedIsolation_->SetDirectory(pDirectory);

        pThirdRho25CorrectedIsolation_=new TH1F( "thirdRho25CorrectedIsolation", "thirdRho25CorrectedIsolation", 120, 0, 2 );
        pThirdRho25CorrectedIsolation_->SetDirectory(pDirectory);

        pDeltaBetaCorrectedIsolationNoZeroing_=new TH1F( "deltaBetaCorrectedIsolationNoZeroing", "deltaBetaCorrectedIsolationNoZeroing", 120, -1, 1 );
        pDeltaBetaCorrectedIsolationNoZeroing_->SetDirectory(pDirectory);

        pRho25CorrectedIsolationNoZeroing_=new TH1F( "rho25CorrectedIsolationNoZeroing", "rho25CorrectedIsolationNoZeroing", 120, -1, 1 );
        pRho25CorrectedIsolationNoZeroing_->SetDirectory(pDirectory);

        pThirdRho25CorrectedIsolationNoZeroing_=new TH1F( "thirdRho25CorrectedIsolationNoZeroing", "thirdRho25CorrectedIsolationNoZeroing", 120, -1, 1 );
        pThirdRho25CorrectedIsolationNoZeroing_->SetDirectory(pDirectory);

        if( createIsolationTree_ )
        {
                pIsolationTree_=new TTree("isolationTree","Data about the lepton isolation");
                pIsolationTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
                pIsolationTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
                pIsolationTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
                pIsolationTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
                pIsolationTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
                pIsolationTree_->Branch("rho25",&rho25_branch_,"rho25/F");
                pIsolationTree_->Branch("pT",&pT_branch_,"pT/F");
                pIsolationTree_->SetDirectory( pDirectory );
        }

        histogramHaveBeenBooked_=true;
}

template<class T_Lepton>
void trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::fill( const T_Lepton& lepton, const VHbbEventAuxInfo* pAuxInfo )
{
        if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::LeptonInfoPlotSet::book() - histograms have not been booked" );

        // Get the additional info about the event if it's available
        int numberOfPrimaryVertices=-1;
        float rho25=-1;
        if( pAuxInfo )
        {
                numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
                rho25=pAuxInfo->puInfo.rho25;
        }

        pEta_->Fill( lepton.p4.Eta() );
        pPhi_->Fill( lepton.p4.Phi() );
        pPt_->Fill( lepton.p4.Pt() );
        pChargedIsolation_->Fill( lepton.pfChaIso );
        pPhotonIsolation_->Fill( lepton.pfPhoIso );
        pNeutralIsolation_->Fill( lepton.pfNeuIso );
        pPileupIsolation_->Fill( lepton.pfChaPUIso );
        pRelativeChargedIsolation_->Fill( lepton.pfChaIso/lepton.p4.Pt() );
        pRelativePhotonIsolation_->Fill( lepton.pfPhoIso/lepton.p4.Pt() );
        pRelativeNeutralIsolation_->Fill( lepton.pfNeuIso/lepton.p4.Pt() );
        pRelativePileupIsolation_->Fill( lepton.pfChaPUIso/lepton.p4.Pt() );
        pRelativeIsolation_->Fill( combinedRelativeIsolation(lepton) );
        pDeltaBetaCorrectedIsolation_->Fill( deltaBetaCorrectedIsolation(lepton) );
        pDeltaBetaCorrectedIsolationNoZeroing_->Fill( deltaBetaCorrectedIsolationNoZeroing(lepton) );
        if( pAuxInfo )
        {
                // Don't bother filling if I couldn't get the rho25 info
                pRho25CorrectedIsolation_->Fill( rho25CorrectedIsolation(lepton,rho25) );
                pRho25CorrectedIsolationNoZeroing_->Fill( rho25CorrectedIsolationNoZeroing(lepton,rho25) );
                pThirdRho25CorrectedIsolation_->Fill( thirdRho25CorrectedIsolation(lepton,rho25) );
                pThirdRho25CorrectedIsolationNoZeroing_->Fill( thirdRho25CorrectedIsolationNoZeroing(lepton,rho25) );
        }

        // If the TTree is not null then isolation data has been requested
        if( pIsolationTree_ )
        {
                numberOfPrimaryVertices_branch_=numberOfPrimaryVertices;
                chargedIsolation_branch_=lepton.pfChaIso;
                photonIsolation_branch_=lepton.pfPhoIso;
                neutralIsolation_branch_=lepton.pfNeuIso;
                pileupIsolation_branch_=lepton.pfChaPUIso;
                rho25_branch_=rho25;
                pT_branch_=lepton.p4.Pt();

                pIsolationTree_->Fill();
        }
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::combinedRelativeIsolation( const T_Lepton& lepton )
{
        return (lepton.pfChaIso+lepton.pfPhoIso+lepton.pfNeuIso)/lepton.p4.Pt();
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::deltaBetaCorrectedIsolation( const T_Lepton& lepton, float deltaBetaFactor )
{
        // deltaBetaFactor defaults to -0.5
        float correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso+deltaBetaFactor*lepton.pfChaPUIso;
        if( correctedNeutralIsolation<0 ) correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso;
        return (lepton.pfChaIso+correctedNeutralIsolation)/lepton.p4.Pt();
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::rho25CorrectedIsolation( const T_Lepton& lepton, float rho25 )
{
        float uncorrectedIsolation=combinedRelativeIsolation(lepton);
        float correctedIsolation=uncorrectedIsolation-(0.4*0.4*M_PI*rho25)/lepton.p4.Pt();
        if( correctedIsolation<0 ) correctedIsolation=0;

        return correctedIsolation;
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::thirdRho25CorrectedIsolation( const T_Lepton& lepton, float rho25 )
{
        float uncorrectedIsolation=combinedRelativeIsolation(lepton);
        float correctedIsolation=uncorrectedIsolation-(0.4*0.4*M_PI*rho25)/3.0/lepton.p4.Pt();
        if( correctedIsolation<0 ) correctedIsolation=0;

        return correctedIsolation;
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::deltaBetaCorrectedIsolationNoZeroing( const T_Lepton& lepton, float deltaBetaFactor )
{
        // deltaBetaFactor defaults to -0.5
        float correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso+deltaBetaFactor*lepton.pfChaPUIso;

        return (lepton.pfChaIso+correctedNeutralIsolation)/lepton.p4.Pt();
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::rho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 )
{
        float uncorrectedIsolation=combinedRelativeIsolation(lepton);
        float correctedIsolation=uncorrectedIsolation-(0.4*0.4*M_PI*rho25)/lepton.p4.Pt();


        return correctedIsolation;
}

template<class T_Lepton>
float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::thirdRho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 )
{
        float uncorrectedIsolation=combinedRelativeIsolation(lepton);
        float correctedIsolation=uncorrectedIsolation-(0.4*0.4*M_PI*rho25)/3.0/lepton.p4.Pt();


        return correctedIsolation;
}

#endif // end of "#ifndef trkupgradeanalysis_LeptonInfoPlotSet_h"
Revision:	1.1
Committed:	Wed Aug 15 22:37:47 2012 UTC (12 years, 8 months ago) by grimes
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	Long overdue commit with several new files
#	User	Rev	Content
1	grimes	1.1	#ifndef trkupgradeanalysis_LeptonInfoPlotSet_h
2			#define trkupgradeanalysis_LeptonInfoPlotSet_h
3
4			#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
5
6			// Forward declarations
7			class TH1F;
8			class TDirectory;
9			class TTree;
10			class VHbbEventAuxInfo;
11
12			namespace trkupgradeanalysis
13			{
14			/** @brief An abstract class that incorporates some of the code used for both electron and muon plotsets.
15			*
16			* Note that this isn't intended to be used on its own, rather as a member for ElectronInfoPlotSet and
17			* MuonInfoPlotSet classes to encapsulate common code between them.
18			*
19			* @author Mark Grimes (mark.grimes@bristol.ac.uk)
20			* @date 02/Jul/2012
21			*/
22			template<class T_Lepton>
23			class LeptonInfoPlotSet
24			{
25			public:
26			/** @brief Only constructor
27			*
28			* createIsolationTree - set whether to record the isolation values in a TTree or not. Default
29			* is off to save space.
30			*/
31			LeptonInfoPlotSet( bool createIsolationTree=false );
32			void book( TDirectory* pDirectory );
33			void fill( const T_Lepton& lepton, const VHbbEventAuxInfo* pAuxInfo=NULL );
34
35			static float combinedRelativeIsolation( const T_Lepton& lepton );
36			static float deltaBetaCorrectedIsolation( const T_Lepton& lepton, float deltaBetaFactor=-0.5 );
37			static float rho25CorrectedIsolation( const T_Lepton& lepton, float rho25 );
38			static float thirdRho25CorrectedIsolation( const T_Lepton& lepton, float rho25 );
39			static float deltaBetaCorrectedIsolationNoZeroing( const T_Lepton& lepton, float deltaBetaFactor=-0.5 );
40			static float rho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 );
41			static float thirdRho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 );
42			private:
43			bool histogramHaveBeenBooked_;
44
45			bool createIsolationTree_; ///< Whether or not to record isolation values in a TTree.
46			TTree* pIsolationTree_; ///< This TTree is only created if createIsolationTree is set to true in the constructor.
47			int numberOfPrimaryVertices_branch_; ///< All these are variables to hold the branch data
48			float chargedIsolation_branch_;
49			float photonIsolation_branch_;
50			float neutralIsolation_branch_;
51			float pileupIsolation_branch_;
52			float rho25_branch_;
53			float pT_branch_;
54
55			TH1F* pEta_;
56			TH1F* pPhi_;
57			TH1F* pPt_;
58			TH1F* pChargedIsolation_;
59			TH1F* pPhotonIsolation_;
60			TH1F* pNeutralIsolation_;
61			TH1F* pPileupIsolation_;
62			TH1F* pRelativeChargedIsolation_;
63			TH1F* pRelativePhotonIsolation_;
64			TH1F* pRelativeNeutralIsolation_;
65			TH1F* pRelativePileupIsolation_;
66			TH1F* pRelativeIsolation_;
67			TH1F* pDeltaBetaCorrectedIsolation_;
68			TH1F* pRho25CorrectedIsolation_;
69			TH1F* pThirdRho25CorrectedIsolation_;
70			TH1F* pDeltaBetaCorrectedIsolationNoZeroing_;
71			TH1F* pRho25CorrectedIsolationNoZeroing_;
72			TH1F* pThirdRho25CorrectedIsolationNoZeroing_;
73			};
74
75			} // end of namespace trkupgradeanalysis
76
77
78
79			//----------------------------------------------------------
80			//
81			// Need the implementation here since it's templated
82			//
83			//----------------------------------------------------------
84			#include <stdexcept>
85			#include <TH1F.h>
86			#include <TTree.h>
87			#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
88
89
90			template<class T_Lepton>
91			trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::LeptonInfoPlotSet( bool createIsolationTree )
92			: histogramHaveBeenBooked_(false), createIsolationTree_(createIsolationTree), pIsolationTree_(NULL)
93			{
94			// No operation besides the initialiser list.
95			}
96
97			template<class T_Lepton>
98			void trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::book( TDirectory* pDirectory )
99			{
100			if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::LeptonInfoPlotSet::book() - histograms have already been booked" );
101
102			//
103			// Note that the root file which TDirectory is part of takes ownership of all
104			// of these objects, so I don't need to (and shouldn't) delete them when I'm
105			// finished.
106			//
107
108
109			pEta_=new TH1F( "eta", "Eta", 60, -3, 3 );
110			pEta_->SetDirectory(pDirectory);
111
112			pPhi_=new TH1F( "phi", "phi", 60, -3.15, 3.15 );
113			pPhi_->SetDirectory(pDirectory);
114
115			pPt_=new TH1F( "pT", "pT", 60, 0, 250 );
116			pPt_->SetDirectory(pDirectory);
117
118			pChargedIsolation_=new TH1F( "chargedIsolation", "chargedIsolation", 60, 0, 30 );
119			pChargedIsolation_->SetDirectory(pDirectory);
120
121			pPhotonIsolation_=new TH1F( "photonIsolation", "photonIsolation", 60, 0, 30 );
122			pPhotonIsolation_->SetDirectory(pDirectory);
123
124			pNeutralIsolation_=new TH1F( "neutralIsolation", "neutralIsolation", 60, 0, 30 );
125			pNeutralIsolation_->SetDirectory(pDirectory);
126
127			pPileupIsolation_=new TH1F( "pileupIsolation", "pileupIsolation", 60, 0, 30 );
128			pPileupIsolation_->SetDirectory(pDirectory);
129
130			pRelativeChargedIsolation_=new TH1F( "relativeChargedIsolation", "relativeChargedIsolation", 60, 0, 1 );
131			pRelativeChargedIsolation_->SetDirectory(pDirectory);
132
133			pRelativePhotonIsolation_=new TH1F( "relativePhotonIsolation", "relativePhotonIsolation", 60, 0, 1 );
134			pRelativePhotonIsolation_->SetDirectory(pDirectory);
135
136			pRelativeNeutralIsolation_=new TH1F( "relativeNeutralIsolation", "relativeNeutralIsolation", 60, 0, 1 );
137			pRelativeNeutralIsolation_->SetDirectory(pDirectory);
138
139			pRelativePileupIsolation_=new TH1F( "relativePileupIsolation", "relativePileupIsolation", 60, 0, 1 );
140			pRelativePileupIsolation_->SetDirectory(pDirectory);
141
142			pRelativeIsolation_=new TH1F( "relativeIsolation", "relativeIsolation", 120, 0, 2 );
143			pRelativeIsolation_->SetDirectory(pDirectory);
144
145			pDeltaBetaCorrectedIsolation_=new TH1F( "deltaBetaCorrectedIsolation", "deltaBetaCorrectedIsolation", 120, 0, 2 );
146			pDeltaBetaCorrectedIsolation_->SetDirectory(pDirectory);
147
148			pRho25CorrectedIsolation_=new TH1F( "rho25CorrectedIsolation", "rho25CorrectedIsolation", 120, 0, 2 );
149			pRho25CorrectedIsolation_->SetDirectory(pDirectory);
150
151			pThirdRho25CorrectedIsolation_=new TH1F( "thirdRho25CorrectedIsolation", "thirdRho25CorrectedIsolation", 120, 0, 2 );
152			pThirdRho25CorrectedIsolation_->SetDirectory(pDirectory);
153
154			pDeltaBetaCorrectedIsolationNoZeroing_=new TH1F( "deltaBetaCorrectedIsolationNoZeroing", "deltaBetaCorrectedIsolationNoZeroing", 120, -1, 1 );
155			pDeltaBetaCorrectedIsolationNoZeroing_->SetDirectory(pDirectory);
156
157			pRho25CorrectedIsolationNoZeroing_=new TH1F( "rho25CorrectedIsolationNoZeroing", "rho25CorrectedIsolationNoZeroing", 120, -1, 1 );
158			pRho25CorrectedIsolationNoZeroing_->SetDirectory(pDirectory);
159
160			pThirdRho25CorrectedIsolationNoZeroing_=new TH1F( "thirdRho25CorrectedIsolationNoZeroing", "thirdRho25CorrectedIsolationNoZeroing", 120, -1, 1 );
161			pThirdRho25CorrectedIsolationNoZeroing_->SetDirectory(pDirectory);
162
163			if( createIsolationTree_ )
164			{
165			pIsolationTree_=new TTree("isolationTree","Data about the lepton isolation");
166			pIsolationTree_->Branch("numberOfPrimaryVertices",&numberOfPrimaryVertices_branch_,"numberOfPrimaryVertices/I");
167			pIsolationTree_->Branch("chargedIsolation",&chargedIsolation_branch_,"chargedIsolation/F");
168			pIsolationTree_->Branch("photonIsolation",&photonIsolation_branch_,"photonIsolation/F");
169			pIsolationTree_->Branch("neutralIsolation",&neutralIsolation_branch_,"neutralIsolation/F");
170			pIsolationTree_->Branch("pileupIsolation",&pileupIsolation_branch_,"pileupIsolation/F");
171			pIsolationTree_->Branch("rho25",&rho25_branch_,"rho25/F");
172			pIsolationTree_->Branch("pT",&pT_branch_,"pT/F");
173			pIsolationTree_->SetDirectory( pDirectory );
174			}
175
176			histogramHaveBeenBooked_=true;
177			}
178
179			template<class T_Lepton>
180			void trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::fill( const T_Lepton& lepton, const VHbbEventAuxInfo* pAuxInfo )
181			{
182			if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::LeptonInfoPlotSet::book() - histograms have not been booked" );
183
184			// Get the additional info about the event if it's available
185			int numberOfPrimaryVertices=-1;
186			float rho25=-1;
187			if( pAuxInfo )
188			{
189			numberOfPrimaryVertices=pAuxInfo->pvInfo.nVertices;
190			rho25=pAuxInfo->puInfo.rho25;
191			}
192
193			pEta_->Fill( lepton.p4.Eta() );
194			pPhi_->Fill( lepton.p4.Phi() );
195			pPt_->Fill( lepton.p4.Pt() );
196			pChargedIsolation_->Fill( lepton.pfChaIso );
197			pPhotonIsolation_->Fill( lepton.pfPhoIso );
198			pNeutralIsolation_->Fill( lepton.pfNeuIso );
199			pPileupIsolation_->Fill( lepton.pfChaPUIso );
200			pRelativeChargedIsolation_->Fill( lepton.pfChaIso/lepton.p4.Pt() );
201			pRelativePhotonIsolation_->Fill( lepton.pfPhoIso/lepton.p4.Pt() );
202			pRelativeNeutralIsolation_->Fill( lepton.pfNeuIso/lepton.p4.Pt() );
203			pRelativePileupIsolation_->Fill( lepton.pfChaPUIso/lepton.p4.Pt() );
204			pRelativeIsolation_->Fill( combinedRelativeIsolation(lepton) );
205			pDeltaBetaCorrectedIsolation_->Fill( deltaBetaCorrectedIsolation(lepton) );
206			pDeltaBetaCorrectedIsolationNoZeroing_->Fill( deltaBetaCorrectedIsolationNoZeroing(lepton) );
207			if( pAuxInfo )
208			{
209			// Don't bother filling if I couldn't get the rho25 info
210			pRho25CorrectedIsolation_->Fill( rho25CorrectedIsolation(lepton,rho25) );
211			pRho25CorrectedIsolationNoZeroing_->Fill( rho25CorrectedIsolationNoZeroing(lepton,rho25) );
212			pThirdRho25CorrectedIsolation_->Fill( thirdRho25CorrectedIsolation(lepton,rho25) );
213			pThirdRho25CorrectedIsolationNoZeroing_->Fill( thirdRho25CorrectedIsolationNoZeroing(lepton,rho25) );
214			}
215
216			// If the TTree is not null then isolation data has been requested
217			if( pIsolationTree_ )
218			{
219			numberOfPrimaryVertices_branch_=numberOfPrimaryVertices;
220			chargedIsolation_branch_=lepton.pfChaIso;
221			photonIsolation_branch_=lepton.pfPhoIso;
222			neutralIsolation_branch_=lepton.pfNeuIso;
223			pileupIsolation_branch_=lepton.pfChaPUIso;
224			rho25_branch_=rho25;
225			pT_branch_=lepton.p4.Pt();
226
227			pIsolationTree_->Fill();
228			}
229			}
230
231			template<class T_Lepton>
232			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::combinedRelativeIsolation( const T_Lepton& lepton )
233			{
234			return (lepton.pfChaIso+lepton.pfPhoIso+lepton.pfNeuIso)/lepton.p4.Pt();
235			}
236
237			template<class T_Lepton>
238			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::deltaBetaCorrectedIsolation( const T_Lepton& lepton, float deltaBetaFactor )
239			{
240			// deltaBetaFactor defaults to -0.5
241			float correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso+deltaBetaFactor*lepton.pfChaPUIso;
242			if( correctedNeutralIsolation<0 ) correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso;
243			return (lepton.pfChaIso+correctedNeutralIsolation)/lepton.p4.Pt();
244			}
245
246			template<class T_Lepton>
247			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::rho25CorrectedIsolation( const T_Lepton& lepton, float rho25 )
248			{
249			float uncorrectedIsolation=combinedRelativeIsolation(lepton);
250			float correctedIsolation=uncorrectedIsolation-(0.40.4M_PI*rho25)/lepton.p4.Pt();
251			if( correctedIsolation<0 ) correctedIsolation=0;
252
253			return correctedIsolation;
254			}
255
256			template<class T_Lepton>
257			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::thirdRho25CorrectedIsolation( const T_Lepton& lepton, float rho25 )
258			{
259			float uncorrectedIsolation=combinedRelativeIsolation(lepton);
260			float correctedIsolation=uncorrectedIsolation-(0.40.4M_PI*rho25)/3.0/lepton.p4.Pt();
261			if( correctedIsolation<0 ) correctedIsolation=0;
262
263			return correctedIsolation;
264			}
265
266			template<class T_Lepton>
267			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::deltaBetaCorrectedIsolationNoZeroing( const T_Lepton& lepton, float deltaBetaFactor )
268			{
269			// deltaBetaFactor defaults to -0.5
270			float correctedNeutralIsolation=lepton.pfPhoIso+lepton.pfNeuIso+deltaBetaFactor*lepton.pfChaPUIso;
271
272			return (lepton.pfChaIso+correctedNeutralIsolation)/lepton.p4.Pt();
273			}
274
275			template<class T_Lepton>
276			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::rho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 )
277			{
278			float uncorrectedIsolation=combinedRelativeIsolation(lepton);
279			float correctedIsolation=uncorrectedIsolation-(0.40.4M_PI*rho25)/lepton.p4.Pt();
280
281
282			return correctedIsolation;
283			}
284
285			template<class T_Lepton>
286			float trkupgradeanalysis::LeptonInfoPlotSet<T_Lepton>::thirdRho25CorrectedIsolationNoZeroing( const T_Lepton& lepton, float rho25 )
287			{
288			float uncorrectedIsolation=combinedRelativeIsolation(lepton);
289			float correctedIsolation=uncorrectedIsolation-(0.40.4M_PI*rho25)/3.0/lepton.p4.Pt();
290
291
292			return correctedIsolation;
293			}
294
295			#endif // end of "#ifndef trkupgradeanalysis_LeptonInfoPlotSet_h"