ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/grimes/VHbbAnalysisCode/src/LightweightAnalyser.cpp

Comparing UserCode/grimes/VHbbAnalysisCode/src/LightweightAnalyser.cpp (file contents):
Revision 1.1 by grimes, Tue Feb 14 01:43:15 2012 UTC vs.
Revision 1.5 by grimes, Wed Aug 15 22:37:47 2012 UTC

#	Line 3 | Line 3
3		#include <sstream>
4
5		#include "TrkUpgradeAnalysis/VHbb/interface/VHbbCandidateCutSets.h"
6	+	#include "TrkUpgradeAnalysis/VHbb/interface/AdditionalJetStudyCutSets.h"
7	+	#include "TrkUpgradeAnalysis/VHbb/interface/tools.h" // required for trkupgradeanalysis::tools::createDirectory
8	+	#include "FWCore/Utilities/interface/EDMException.h"
9
10		#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
11		#include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc" // Not entirely sure why this is included and not linked to
#	Line 17 | Line 20
20
21
22		const bool trkupgradeanalysis::LightweightAnalyser::verbose_=false;
23	<	const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdZ_=20;
24	<	const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdW_=20;
23	>	const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdZ_=30;
24	>	const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdW_=30;
25		const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsZ_=false;
26		const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsW_=true;
27
25	–	namespace // Use the unnamed namespace
26	–	{
27	–	TDirectory* createDirectory( const std::string& fullPath, TDirectory* pParent )
28	–	{
29	–	if( pParent==NULL ) throw std::runtime_error( "The parent directory is a Null pointer" );
30	–
31	–	TDirectory* pSubDirectory=pParent;
32	–	size_t currentPosition=0;
33	–	size_t nextSlash;
34	–	do
35	–	{
36	–	nextSlash=fullPath.find_first_of('/', currentPosition );
37	–	std::string directoryName=fullPath.substr(currentPosition,nextSlash-currentPosition);
38	–	currentPosition=nextSlash+1;
39	–
40	–	TDirectory* pNextSubDirectory=pSubDirectory->GetDirectory( directoryName.c_str() );
41	–	if( pNextSubDirectory==NULL ) pNextSubDirectory=pSubDirectory->mkdir( directoryName.c_str() );
42	–	if( pNextSubDirectory==NULL ) throw std::runtime_error( "Couldn't create the root directory \""+directoryName+"\"" );
43	–	pSubDirectory=pNextSubDirectory;
44	–
45	–	} while( nextSlash!=std::string::npos );
46	–
47	–	return pSubDirectory;
48	–	}
49	–	}
28
29		trkupgradeanalysis::LightweightAnalyser::LightweightAnalyser( const std::string& outputFilename )
30		: zFinder_( verbose_, jetPtThresholdZ_, useHighestPtHiggsZ_ ),
#	Line 58 | Line 36 | trkupgradeanalysis::LightweightAnalyser:
36		// All of these are empty smart pointers, so I need to call reset on them
37		pOutputFile_.reset( new TFile( outputFilename.c_str(), "RECREATE" ) );
38
39	<	cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::SignalSelectionZee(120)) );
40	<	cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::SignalSelectionZmumu(120)) );
41	<	cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::SignalSelectionWen(120)) );
42	<	cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::SignalSelectionWmun(120)) );
43	<	cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::SignalSelectionZmumuWithoutAdditionalJetsCut(120)) );
44	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VlightRegionHWmun) );
45	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VlightRegionHWen) );
46	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VlightRegionHZmumu) );
47	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VlightRegionHZee) );
70	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::TTbarRegionHWmun) );
71	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::TTbarRegionHWen) );
72	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::TTbarRegionHZmumu) );
73	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::TTbarRegionHZee) );
74	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VbbRegionHWmun) );
75	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VbbRegionHWen) );
76	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VbbRegionHZmumu) );
77	<	//      cutsToApply_.push_back( boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet>(new trkupgradeanalysis::VbbRegionHZee) );
39	>
40	>	using namespace trkupgradeanalysis; // Don't really need this but I like to be explicit (ooh err..)
41	>	//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZee(120)) );
42	>	//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(120)) );
43	>	cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new AllVariables) );
44	>	//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionStdGeom50PUZmumuWithoutAdditionalJetsCut(125)) );
45	>	//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionPhase150PUZmumuWithoutAdditionalJetsCut(125)) );
46	>	cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new InitialSignalSelectionZee) );
47	>	cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new InitialSignalSelectionZmumu) );
48		}
49
50		trkupgradeanalysis::LightweightAnalyser::~LightweightAnalyser()
#	Line 108 | Line 78 | void trkupgradeanalysis::LightweightAnal
78		void trkupgradeanalysis::LightweightAnalyser::bookHistograms( trkupgradeanalysis::LightweightAnalyser::AllPlots& plots, const std::string& directoryName )
79		{
80		TDirectory* pRootDirectory=pOutputFile_->GetDirectory("/");
81	<	TDirectory* pDirectory=createDirectory( directoryName, pRootDirectory );
81	>	TDirectory* pDirectory=trkupgradeanalysis::tools::createDirectory( directoryName, pRootDirectory );
82		//      TDirectory* pDirectory=pRootDirectory->GetDirectory( directoryName.c_str() );
83		//      if( pDirectory==NULL ) pDirectory=pRootDirectory->mkdir( directoryName.c_str() );
84
#	Line 126 | Line 96 | void trkupgradeanalysis::LightweightAnal
96		pSubDirectory=pDirectory->mkdir( "ElectronInfos" );
97		plots.allElectrons_.book( pSubDirectory );
98
99	+	pSubDirectory=pDirectory->mkdir( "HiggsCandidateJets" );
100	+	plots.higgsCandidateJets_.book( pSubDirectory );
101	+
102	+	pSubDirectory=pDirectory->mkdir( "AdditionalJets" );
103	+	plots.additionalJets_.book( pSubDirectory );
104	+
105	+	pSubDirectory=pDirectory->mkdir( "SimpleJets2" );
106	+	plots.simpleJets2_.book( pSubDirectory );
107	+
108	+	pSubDirectory=pDirectory->mkdir( "IsolationStudy" );
109	+	plots.isolationStudy_.book( pSubDirectory );
110	+
111		pSubDirectory=pDirectory->mkdir( "MonteCarloTruth" );
112		plots.monteCarloInfo_.book( pSubDirectory );
113
#	Line 142 | Line 124 | void trkupgradeanalysis::LightweightAnal
124		}
125		}
126
127	<	void trkupgradeanalysis::LightweightAnalyser::fillAllPlotsStructure( AllPlots* pPlotsToFill, const VHbbEvent& vhbbEvent, const std::vector<VHbbCandidate>& zCandidates, const std::vector<VHbbCandidate>& wCandidates )
127	>	void trkupgradeanalysis::LightweightAnalyser::fillAllPlotsStructure( AllPlots* pPlotsToFill, const VHbbEvent& vhbbEvent, const VHbbEventAuxInfo* pVHbbEventAuxInfo, const std::vector<VHbbCandidate>& zCandidates, const std::vector<VHbbCandidate>& wCandidates )
128		{
129	<	pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo );
129	>	pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo, pVHbbEventAuxInfo );
130		pPlotsToFill->allElectrons_.fill( vhbbEvent.eleInfo );
131	+	if( pVHbbEventAuxInfo ) pPlotsToFill->monteCarloInfo_.fill( *pVHbbEventAuxInfo );
132	+	pPlotsToFill->isolationStudy_.fill( vhbbEvent, pVHbbEventAuxInfo ); // This class checks if is null itself
133	+	pPlotsToFill->simpleJets2_.fill( vhbbEvent.simpleJets2, pVHbbEventAuxInfo );
134
135		for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
136		{
137		pPlotsToFill->candidateHistogramsForAllEvents_.fill( *iCandidate );
138	+	pPlotsToFill->higgsCandidateJets_.fill( iCandidate->H.jets, pVHbbEventAuxInfo );
139	+	pPlotsToFill->additionalJets_.fill( iCandidate->additionalJets, pVHbbEventAuxInfo );
140
141		// Loop over all of the cuts. If this candidate passes the cuts then fill the plots
142		for( std::vector<CutSetPlotSet>::iterator iCutPlotSet=pPlotsToFill->cutCollectionPlotSets_.begin();
#	Line 180 | Line 167 | void trkupgradeanalysis::LightweightAnal
167		fwlite::Handle<VHbbEventAuxInfo> vhbbAuxInfoHandle;
168		vhbbAuxInfoHandle.getByLabel( event, "HbbAnalyzerNew" );
169		const VHbbEventAuxInfo* pVHbbEventAuxInfo=vhbbAuxInfoHandle.product();
170	<	if( pVHbbEventAuxInfo == 0 )
184	<	{
185	<	std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
186	<	return;
187	<	}
170	>	if( pVHbbEventAuxInfo == 0 ) std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
171		else
172		{
190	–	pPlotsForEverything_->allEventTypes_.monteCarloInfo_.fill( *pVHbbEventAuxInfo );
173		std::string eventType=eventTypeFromMC( *pVHbbEventAuxInfo );
174
175		std::map<std::string,AllPlots>::iterator iPlotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_.find(eventType);
#	Line 204 | Line 186 | void trkupgradeanalysis::LightweightAnal
186
187		++numberOfEvents_;
188
189	+	bool applyMuonIsolation=false;
190	+	// Apply isolation to the muons. Need to modify so use a very naughty const cast
191	+	if( applyMuonIsolation )
192	+	{
193	+	std::vector<VHbbEvent::MuonInfo>& muonInfos=const_cast< std::vector<VHbbEvent::MuonInfo>& >( pVHbbEvent->muInfo );
194	+	bool allMuonsPassIsolationCut=false;
195	+	while( !allMuonsPassIsolationCut )
196	+	{
197	+	// I'm worried about looping using invalidated iterators, so whenever I erase an element I'll
198	+	// start looping from the start. Not the most efficient way but these vectors typically have
199	+	// at most 3 entries.
200	+	std::vector<VHbbEvent::MuonInfo>::iterator iMuon;
201	+	for( iMuon=muonInfos.begin(); iMuon!=muonInfos.end(); ++iMuon )
202	+	{
203	+	float isolation=MuonInfoPlotSet::thirdRho25CorrectedIsolation( *iMuon, pVHbbEventAuxInfo->puInfo.rho25 );
204	+	if( isolation>=0.15 ) break;
205	+	}
206	+	if( iMuon!=muonInfos.end() ) muonInfos.erase( iMuon );
207	+	else allMuonsPassIsolationCut=true;
208	+	}
209	+	}
210	+
211		zFinder_.run( pVHbbEvent, zCandidates );
212		wFinder_.run( pVHbbEvent, wCandidates );
213
214	<	//      fwlite::Handle<edm::View<pat::Muon> > muonHandle;
215	<	//      muonHandle.getByLabel( event, "muons" );
216	<	//      const edm::View<pat::Muon>& muons=*muonHandle;
213	<	//      std::cout << "muons.size()=" << muons.size() << std::endl;
214	<
215	<	fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, zCandidates, wCandidates );
216	<	if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, zCandidates, wCandidates );
217	<
218	<	//      // This loops over all of the muons in the event and plots some things about them
219	<	//      pCurrentPlots_->allMuons_.fill( pVHbbEvent->muInfo );
220	<	//      pCurrentPlots_->allElectrons_.fill( pVHbbEvent->eleInfo );
221	<	//
222	<	//      for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
223	<	//      {
224	<	//              pCurrentPlots_->candidateHistogramsForAllEvents_.fill( *iCandidate );
225	<	//
226	<	//              // Loop over all of the cuts. If this candidate passes the cuts then fill the plots
227	<	//              for( std::vector<CutSetPlotSet>::iterator iCutPlotSet=pCurrentPlots_->cutCollectionPlotSets_.begin();
228	<	//                              iCutPlotSet!=pCurrentPlots_->cutCollectionPlotSets_.end(); ++iCutPlotSet )
229	<	//              {
230	<	//                      iCutPlotSet->fill( *iCandidate );
231	<	//              }
232	<	//      }
214	>	fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
215	>	if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
216	>
217
218		if( !zCandidates.empty() || !wCandidates.empty() ) ++numberOfEventsWithAtLeastOneZOrWCandidate_;
219
236	–	//      pCurrentPlots_->pNumberOfCandidates_->Fill( zCandidates.size(), wCandidates.size() );
220		}
221
222	<	void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile )
222	>	void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile, size_t numberOfEvents )
223		{
224		if( pPlotsForEverything_==NULL ) changeDirectory("allEvents");
225
226		fwlite::Event event( pInputFile );
227	+
228	+	size_t eventCount=0;
229		for( event.toBegin(); !event.atEnd(); ++event )
230		{
231	<	processEvent( event );
231	>	try
232	>	{
233	>	processEvent( event );
234	>	++eventCount;
235	>	}
236	>	catch( edm::Exception& error )
237	>	{
238	>	std::cerr << error.what() << std::endl;
239	>	}
240	>
241	>	if( numberOfEvents!=0 && eventCount>=numberOfEvents )
242	>	{
243	>	std::cout << "Breaking after " << eventCount << " events" << std::endl;
244	>	break;
245	>	}
246		}
247
248		std::cout << "There were " << numberOfEventsWithAtLeastOneZOrWCandidate_ << " events with a Z or W candidate out of " << numberOfEvents_ << " events." << std::endl;
#	Line 286 | Line 285 | std::string trkupgradeanalysis::Lightwei
285		eventStringStream << ") ";
286		}
287
288	<	return eventStringStream.str();
288	>	if( eventStringStream.str().empty() ) return "otherBackground";
289	>	else return eventStringStream.str();
290		}
291

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines