VHbbAnalysisCode/src/LightweightAnalyser.cpp

#include "TrkUpgradeAnalysis/VHbb/interface/LightweightAnalyser.h"

#include <sstream>

#include "TrkUpgradeAnalysis/VHbb/interface/VHbbCandidateCutSets.h"
#include "TrkUpgradeAnalysis/VHbb/interface/AdditionalJetStudyCutSets.h"
#include "TrkUpgradeAnalysis/VHbb/interface/tools.h" // required for trkupgradeanalysis::tools::createDirectory
#include "FWCore/Utilities/interface/EDMException.h"

#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
#include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc" // Not entirely sure why this is included and not linked to
#include "DataFormats/FWLite/interface/Event.h"
#include "DataFormats/FWLite/interface/Handle.h"
#include "DataFormats/PatCandidates/interface/Muon.h"


#include <TFile.h>
#include <TH1.h>
#include <TH2.h>


const bool trkupgradeanalysis::LightweightAnalyser::verbose_=false;
const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdZ_=30;
const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdW_=30;
const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsZ_=false;
const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsW_=true;


trkupgradeanalysis::LightweightAnalyser::LightweightAnalyser( const std::string& outputFilename )
        : zFinder_( verbose_, jetPtThresholdZ_, useHighestPtHiggsZ_ ),
          wFinder_( verbose_, jetPtThresholdW_, useHighestPtHiggsW_ ),
          numberOfEvents_(0),
          numberOfEventsWithAtLeastOneZOrWCandidate_(0),
          pPlotsForEverything_(NULL)
{
        // All of these are empty smart pointers, so I need to call reset on them
        pOutputFile_.reset( new TFile( outputFilename.c_str(), "RECREATE" ) );


        using namespace trkupgradeanalysis; // Don't really need this but I like to be explicit (ooh err..)
//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZee(120)) );
//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(120)) );
        cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new AllVariables) );
//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionStdGeom50PUZmumuWithoutAdditionalJetsCut(125)) );
//      cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionPhase150PUZmumuWithoutAdditionalJetsCut(125)) );
        cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new InitialSignalSelectionZee) );
        cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new InitialSignalSelectionZmumu) );
}

trkupgradeanalysis::LightweightAnalyser::~LightweightAnalyser()
{
        save();
        pOutputFile_->Close();
}

void trkupgradeanalysis::LightweightAnalyser::changeDirectory( const std::string& directory )
{
        std::map<std::string,PlotsSplitByMCType>::iterator iPlots=plotsForEventTypes_.find(directory);

        if( iPlots==plotsForEventTypes_.end() )
        {
                pPlotsForEverything_=&plotsForEventTypes_[directory]; // This will create a new entry in the map

                // Only book the histograms for all MC event types now. I don't know what MC event types there
                // until I look at the events so I need to wait until then
                bookHistograms( pPlotsForEverything_->allEventTypes_, directory+"/allEvents" );
        }
        else pPlotsForEverything_=&iPlots->second; // This directory has already been created and the histograms booked, so reuse those

        currentDirectory_=directory;
}

void trkupgradeanalysis::LightweightAnalyser::save()
{
        pOutputFile_->Write(0,TObject::kOverwrite);
}

void trkupgradeanalysis::LightweightAnalyser::bookHistograms( trkupgradeanalysis::LightweightAnalyser::AllPlots& plots, const std::string& directoryName )
{
        TDirectory* pRootDirectory=pOutputFile_->GetDirectory("/");
        TDirectory* pDirectory=trkupgradeanalysis::tools::createDirectory( directoryName, pRootDirectory );
//      TDirectory* pDirectory=pRootDirectory->GetDirectory( directoryName.c_str() );
//      if( pDirectory==NULL ) pDirectory=pRootDirectory->mkdir( directoryName.c_str() );

        TDirectory* pSubDirectory=pDirectory->mkdir( "VHbbCandidates" );

        plots.pNumberOfCandidates_=new TH2F( "numberOfCandidates","Number of candidates",5,-0.5,4.5,5,-0.5,4.5);
        plots.pNumberOfCandidates_->SetDirectory(pSubDirectory);
        plots.candidateHistogramsForAllEvents_.book( pSubDirectory );
        plots.pNumberOfCandidates_->GetXaxis()->SetTitle( "Z candidates" );
        plots.pNumberOfCandidates_->GetYaxis()->SetTitle( "W candidates" );

        pSubDirectory=pDirectory->mkdir( "MuonInfos" );
        plots.allMuons_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "ElectronInfos" );
        plots.allElectrons_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "HiggsCandidateJets" );
        plots.higgsCandidateJets_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "AdditionalJets" );
        plots.additionalJets_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "SimpleJets2" );
        plots.simpleJets2_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "IsolationStudy" );
        plots.isolationStudy_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "MonteCarloTruth" );
        plots.monteCarloInfo_.book( pSubDirectory );

        pSubDirectory=pDirectory->mkdir( "Cuts" );

        for( std::vector< boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet> >::const_iterator iCutSet=cutsToApply_.begin();
                        iCutSet!=cutsToApply_.end(); ++iCutSet )
        {
                TDirectory* pSubSubDirectory=pSubDirectory->mkdir( (*iCutSet)->name().c_str() );

                CutSetPlotSet newPlotSet( *iCutSet );
                newPlotSet.book( pSubSubDirectory );
                plots.cutCollectionPlotSets_.push_back( newPlotSet );
        }
}

void trkupgradeanalysis::LightweightAnalyser::fillAllPlotsStructure( AllPlots* pPlotsToFill, const VHbbEvent& vhbbEvent, const VHbbEventAuxInfo* pVHbbEventAuxInfo, const std::vector<VHbbCandidate>& zCandidates, const std::vector<VHbbCandidate>& wCandidates )
{
        pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo, pVHbbEventAuxInfo );
        pPlotsToFill->allElectrons_.fill( vhbbEvent.eleInfo );
        if( pVHbbEventAuxInfo ) pPlotsToFill->monteCarloInfo_.fill( *pVHbbEventAuxInfo );
        pPlotsToFill->isolationStudy_.fill( vhbbEvent, pVHbbEventAuxInfo ); // This class checks if is null itself
        pPlotsToFill->simpleJets2_.fill( vhbbEvent.simpleJets2, pVHbbEventAuxInfo );

        for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
        {
                pPlotsToFill->candidateHistogramsForAllEvents_.fill( *iCandidate );
                pPlotsToFill->higgsCandidateJets_.fill( iCandidate->H.jets, pVHbbEventAuxInfo );
                pPlotsToFill->additionalJets_.fill( iCandidate->additionalJets, pVHbbEventAuxInfo );

                // Loop over all of the cuts. If this candidate passes the cuts then fill the plots
                for( std::vector<CutSetPlotSet>::iterator iCutPlotSet=pPlotsToFill->cutCollectionPlotSets_.begin();
                                iCutPlotSet!=pPlotsToFill->cutCollectionPlotSets_.end(); ++iCutPlotSet )
                {
                        iCutPlotSet->fill( *iCandidate );
                }
        }

        pPlotsToFill->pNumberOfCandidates_->Fill( zCandidates.size(), wCandidates.size() );
}

void trkupgradeanalysis::LightweightAnalyser::processEvent( const fwlite::Event& event )
{
        std::vector<VHbbCandidate> zCandidates;
        std::vector<VHbbCandidate> wCandidates;

        fwlite::Handle<VHbbEvent> vhbbHandle;
        vhbbHandle.getByLabel( event, "HbbAnalyzerNew" );
        const VHbbEvent* pVHbbEvent=vhbbHandle.product();
        if( pVHbbEvent == 0 )
        {
                std::cerr << "Couldn't get the VHbbEvent" << std::endl;
                return;
        }

        AllPlots* pCurrentPlotsForThisEventType=NULL;
        fwlite::Handle<VHbbEventAuxInfo> vhbbAuxInfoHandle;
        vhbbAuxInfoHandle.getByLabel( event, "HbbAnalyzerNew" );
        const VHbbEventAuxInfo* pVHbbEventAuxInfo=vhbbAuxInfoHandle.product();
        if( pVHbbEventAuxInfo == 0 ) std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
        else
        {
                std::string eventType=eventTypeFromMC( *pVHbbEventAuxInfo );

                std::map<std::string,AllPlots>::iterator iPlotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_.find(eventType);

                if( iPlotsForThisEventType==pPlotsForEverything_->mcEventTypePlots_.end() )
                {
                        // This event type hasn't been encountered before so I need to book the histograms
                        AllPlots& plotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_[eventType];
                        bookHistograms( plotsForThisEventType, currentDirectory_+"/"+eventType );
                        pCurrentPlotsForThisEventType=&plotsForThisEventType;
                }
                else pCurrentPlotsForThisEventType=&(iPlotsForThisEventType->second);
        }

        ++numberOfEvents_;

        bool applyMuonIsolation=false;
        // Apply isolation to the muons. Need to modify so use a very naughty const cast
        if( applyMuonIsolation )
        {
                std::vector<VHbbEvent::MuonInfo>& muonInfos=const_cast< std::vector<VHbbEvent::MuonInfo>& >( pVHbbEvent->muInfo );
                bool allMuonsPassIsolationCut=false;
                while( !allMuonsPassIsolationCut )
                {
                        // I'm worried about looping using invalidated iterators, so whenever I erase an element I'll
                        // start looping from the start. Not the most efficient way but these vectors typically have
                        // at most 3 entries.
                        std::vector<VHbbEvent::MuonInfo>::iterator iMuon;
                        for( iMuon=muonInfos.begin(); iMuon!=muonInfos.end(); ++iMuon )
                        {
                                float isolation=MuonInfoPlotSet::thirdRho25CorrectedIsolation( *iMuon, pVHbbEventAuxInfo->puInfo.rho25 );
                                if( isolation>=0.15 ) break;
                        }
                        if( iMuon!=muonInfos.end() ) muonInfos.erase( iMuon );
                        else allMuonsPassIsolationCut=true;
                }
        }

        zFinder_.run( pVHbbEvent, zCandidates );
        wFinder_.run( pVHbbEvent, wCandidates );

        fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
        if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );


        if( !zCandidates.empty() || !wCandidates.empty() ) ++numberOfEventsWithAtLeastOneZOrWCandidate_;

}

void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile, size_t numberOfEvents )
{
        if( pPlotsForEverything_==NULL ) changeDirectory("allEvents");

        fwlite::Event event( pInputFile );

        size_t eventCount=0;
        for( event.toBegin(); !event.atEnd(); ++event )
        {
                try
                {
                        processEvent( event );
                        ++eventCount;
                }
                catch( edm::Exception& error )
                {
                        std::cerr << error.what() << std::endl;
                }

                if( numberOfEvents!=0 && eventCount>=numberOfEvents )
                {
                        std::cout << "Breaking after " << eventCount << " events" << std::endl;
                        break;
                }
        }

        std::cout << "There were " << numberOfEventsWithAtLeastOneZOrWCandidate_ << " events with a Z or W candidate out of " << numberOfEvents_ << " events." << std::endl;
}

std::string trkupgradeanalysis::LightweightAnalyser::eventTypeFromMC( const VHbbEventAuxInfo& eventAuxInfo )
{
        std::stringstream eventStringStream;

        std::vector<VHbbEventAuxInfo::ParticleMCInfo> zMCInfo=eventAuxInfo.mcZ;
        std::vector<VHbbEventAuxInfo::ParticleMCInfo> higgsMCInfo=eventAuxInfo.mcH;

        for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator higgsInfo=higgsMCInfo.begin(); higgsInfo!=higgsMCInfo.end(); ++higgsInfo )
        {
                if( higgsInfo->dauid.size()==1 && std::abs(higgsInfo->dauid[0])==25 ) continue;

                eventStringStream << "H->(";
                for( std::vector<int>::const_iterator daughterID=higgsInfo->dauid.begin(); daughterID!=higgsInfo->dauid.end(); ++daughterID )
                {
                        if( daughterID!=higgsInfo->dauid.begin() ) eventStringStream << ",";
                        eventStringStream << *daughterID;
                }
                eventStringStream << ") ";
        }

        for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator zInfo=zMCInfo.begin(); zInfo!=zMCInfo.end(); ++zInfo )
        {
                if( zInfo->dauid.empty() ) continue;
                if( zInfo->dauid.size()==1 && std::abs(zInfo->dauid[0])==23 ) continue;
                if( zInfo->dauid.size()==2 && std::abs(zInfo->dauid[0])==23 && std::abs(zInfo->dauid[1])==25 ) continue;

                eventStringStream << "Z->(";
                for( std::vector<int>::const_iterator daughterID=zInfo->dauid.begin(); daughterID!=zInfo->dauid.end(); ++daughterID )
                {
                        if( *daughterID==23 ) continue;
                        if( zInfo->dauid.size()>2 && *daughterID==22 ) continue;
                        if( daughterID!=zInfo->dauid.begin() ) eventStringStream << ",";
                        eventStringStream << *daughterID;
                }
                eventStringStream << ") ";
        }

        if( eventStringStream.str().empty() ) return "otherBackground";
        else return eventStringStream.str();
}

Revision:	1.5
Committed:	Wed Aug 15 22:37:47 2012 UTC (12 years, 8 months ago) by grimes
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.4:	+65 -73 lines
Log Message:	Long overdue commit with several new files
#	User	Rev	Content
1	grimes	1.1	#include "TrkUpgradeAnalysis/VHbb/interface/LightweightAnalyser.h"
2
3			#include <sstream>
4
5			#include "TrkUpgradeAnalysis/VHbb/interface/VHbbCandidateCutSets.h"
6	grimes	1.2	#include "TrkUpgradeAnalysis/VHbb/interface/AdditionalJetStudyCutSets.h"
7	grimes	1.5	#include "TrkUpgradeAnalysis/VHbb/interface/tools.h" // required for trkupgradeanalysis::tools::createDirectory
8			#include "FWCore/Utilities/interface/EDMException.h"
9	grimes	1.1
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
12			#include "DataFormats/FWLite/interface/Event.h"
13			#include "DataFormats/FWLite/interface/Handle.h"
14			#include "DataFormats/PatCandidates/interface/Muon.h"
15
16
17			#include <TFile.h>
18			#include <TH1.h>
19			#include <TH2.h>
20
21
22			const bool trkupgradeanalysis::LightweightAnalyser::verbose_=false;
23	grimes	1.5	const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdZ_=30;
24			const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdW_=30;
25	grimes	1.1	const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsZ_=false;
26			const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsW_=true;
27
28
29			trkupgradeanalysis::LightweightAnalyser::LightweightAnalyser( const std::string& outputFilename )
30			: zFinder_( verbose_, jetPtThresholdZ_, useHighestPtHiggsZ_ ),
31			wFinder_( verbose_, jetPtThresholdW_, useHighestPtHiggsW_ ),
32			numberOfEvents_(0),
33			numberOfEventsWithAtLeastOneZOrWCandidate_(0),
34			pPlotsForEverything_(NULL)
35			{
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	grimes	1.2
40			using namespace trkupgradeanalysis; // Don't really need this but I like to be explicit (ooh err..)
41	grimes	1.5	// 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	grimes	1.1	}
49
50			trkupgradeanalysis::LightweightAnalyser::~LightweightAnalyser()
51			{
52			save();
53			pOutputFile_->Close();
54			}
55
56			void trkupgradeanalysis::LightweightAnalyser::changeDirectory( const std::string& directory )
57			{
58			std::map<std::string,PlotsSplitByMCType>::iterator iPlots=plotsForEventTypes_.find(directory);
59
60			if( iPlots==plotsForEventTypes_.end() )
61			{
62			pPlotsForEverything_=&plotsForEventTypes_[directory]; // This will create a new entry in the map
63
64			// Only book the histograms for all MC event types now. I don't know what MC event types there
65			// until I look at the events so I need to wait until then
66			bookHistograms( pPlotsForEverything_->allEventTypes_, directory+"/allEvents" );
67			}
68			else pPlotsForEverything_=&iPlots->second; // This directory has already been created and the histograms booked, so reuse those
69
70			currentDirectory_=directory;
71			}
72
73			void trkupgradeanalysis::LightweightAnalyser::save()
74			{
75			pOutputFile_->Write(0,TObject::kOverwrite);
76			}
77
78			void trkupgradeanalysis::LightweightAnalyser::bookHistograms( trkupgradeanalysis::LightweightAnalyser::AllPlots& plots, const std::string& directoryName )
79			{
80			TDirectory* pRootDirectory=pOutputFile_->GetDirectory("/");
81	grimes	1.5	TDirectory* pDirectory=trkupgradeanalysis::tools::createDirectory( directoryName, pRootDirectory );
82	grimes	1.1	// TDirectory* pDirectory=pRootDirectory->GetDirectory( directoryName.c_str() );
83			// if( pDirectory==NULL ) pDirectory=pRootDirectory->mkdir( directoryName.c_str() );
84
85			TDirectory* pSubDirectory=pDirectory->mkdir( "VHbbCandidates" );
86
87			plots.pNumberOfCandidates_=new TH2F( "numberOfCandidates","Number of candidates",5,-0.5,4.5,5,-0.5,4.5);
88			plots.pNumberOfCandidates_->SetDirectory(pSubDirectory);
89			plots.candidateHistogramsForAllEvents_.book( pSubDirectory );
90			plots.pNumberOfCandidates_->GetXaxis()->SetTitle( "Z candidates" );
91			plots.pNumberOfCandidates_->GetYaxis()->SetTitle( "W candidates" );
92
93			pSubDirectory=pDirectory->mkdir( "MuonInfos" );
94			plots.allMuons_.book( pSubDirectory );
95
96			pSubDirectory=pDirectory->mkdir( "ElectronInfos" );
97			plots.allElectrons_.book( pSubDirectory );
98
99	grimes	1.5	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	grimes	1.4	pSubDirectory=pDirectory->mkdir( "IsolationStudy" );
109			plots.isolationStudy_.book( pSubDirectory );
110
111	grimes	1.1	pSubDirectory=pDirectory->mkdir( "MonteCarloTruth" );
112			plots.monteCarloInfo_.book( pSubDirectory );
113
114			pSubDirectory=pDirectory->mkdir( "Cuts" );
115
116			for( std::vector< boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet> >::const_iterator iCutSet=cutsToApply_.begin();
117			iCutSet!=cutsToApply_.end(); ++iCutSet )
118			{
119			TDirectory* pSubSubDirectory=pSubDirectory->mkdir( (*iCutSet)->name().c_str() );
120
121			CutSetPlotSet newPlotSet( *iCutSet );
122			newPlotSet.book( pSubSubDirectory );
123			plots.cutCollectionPlotSets_.push_back( newPlotSet );
124			}
125			}
126
127	grimes	1.4	void trkupgradeanalysis::LightweightAnalyser::fillAllPlotsStructure( AllPlots* pPlotsToFill, const VHbbEvent& vhbbEvent, const VHbbEventAuxInfo* pVHbbEventAuxInfo, const std::vector<VHbbCandidate>& zCandidates, const std::vector<VHbbCandidate>& wCandidates )
128	grimes	1.1	{
129	grimes	1.5	pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo, pVHbbEventAuxInfo );
130	grimes	1.1	pPlotsToFill->allElectrons_.fill( vhbbEvent.eleInfo );
131	grimes	1.4	if( pVHbbEventAuxInfo ) pPlotsToFill->monteCarloInfo_.fill( *pVHbbEventAuxInfo );
132			pPlotsToFill->isolationStudy_.fill( vhbbEvent, pVHbbEventAuxInfo ); // This class checks if is null itself
133	grimes	1.5	pPlotsToFill->simpleJets2_.fill( vhbbEvent.simpleJets2, pVHbbEventAuxInfo );
134	grimes	1.1
135			for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
136			{
137			pPlotsToFill->candidateHistogramsForAllEvents_.fill( *iCandidate );
138	grimes	1.5	pPlotsToFill->higgsCandidateJets_.fill( iCandidate->H.jets, pVHbbEventAuxInfo );
139			pPlotsToFill->additionalJets_.fill( iCandidate->additionalJets, pVHbbEventAuxInfo );
140	grimes	1.1
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();
143			iCutPlotSet!=pPlotsToFill->cutCollectionPlotSets_.end(); ++iCutPlotSet )
144			{
145			iCutPlotSet->fill( *iCandidate );
146			}
147			}
148
149			pPlotsToFill->pNumberOfCandidates_->Fill( zCandidates.size(), wCandidates.size() );
150			}
151
152			void trkupgradeanalysis::LightweightAnalyser::processEvent( const fwlite::Event& event )
153			{
154			std::vector<VHbbCandidate> zCandidates;
155			std::vector<VHbbCandidate> wCandidates;
156
157			fwlite::Handle<VHbbEvent> vhbbHandle;
158			vhbbHandle.getByLabel( event, "HbbAnalyzerNew" );
159			const VHbbEvent* pVHbbEvent=vhbbHandle.product();
160			if( pVHbbEvent == 0 )
161			{
162			std::cerr << "Couldn't get the VHbbEvent" << std::endl;
163			return;
164			}
165
166			AllPlots* pCurrentPlotsForThisEventType=NULL;
167			fwlite::Handle<VHbbEventAuxInfo> vhbbAuxInfoHandle;
168			vhbbAuxInfoHandle.getByLabel( event, "HbbAnalyzerNew" );
169			const VHbbEventAuxInfo* pVHbbEventAuxInfo=vhbbAuxInfoHandle.product();
170	grimes	1.2	if( pVHbbEventAuxInfo == 0 ) std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
171	grimes	1.1	else
172			{
173			std::string eventType=eventTypeFromMC( *pVHbbEventAuxInfo );
174
175			std::map<std::string,AllPlots>::iterator iPlotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_.find(eventType);
176
177			if( iPlotsForThisEventType==pPlotsForEverything_->mcEventTypePlots_.end() )
178			{
179			// This event type hasn't been encountered before so I need to book the histograms
180			AllPlots& plotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_[eventType];
181			bookHistograms( plotsForThisEventType, currentDirectory_+"/"+eventType );
182			pCurrentPlotsForThisEventType=&plotsForThisEventType;
183			}
184			else pCurrentPlotsForThisEventType=&(iPlotsForThisEventType->second);
185			}
186
187			++numberOfEvents_;
188
189	grimes	1.5	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	grimes	1.1	zFinder_.run( pVHbbEvent, zCandidates );
212			wFinder_.run( pVHbbEvent, wCandidates );
213
214	grimes	1.4	fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
215			if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
216	grimes	1.1
217
218			if( !zCandidates.empty() \|\| !wCandidates.empty() ) ++numberOfEventsWithAtLeastOneZOrWCandidate_;
219
220			}
221
222	grimes	1.5	void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile, size_t numberOfEvents )
223	grimes	1.1	{
224			if( pPlotsForEverything_==NULL ) changeDirectory("allEvents");
225
226			fwlite::Event event( pInputFile );
227	grimes	1.5
228			size_t eventCount=0;
229	grimes	1.1	for( event.toBegin(); !event.atEnd(); ++event )
230			{
231	grimes	1.5	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	grimes	1.1	}
247
248			std::cout << "There were " << numberOfEventsWithAtLeastOneZOrWCandidate_ << " events with a Z or W candidate out of " << numberOfEvents_ << " events." << std::endl;
249			}
250
251			std::string trkupgradeanalysis::LightweightAnalyser::eventTypeFromMC( const VHbbEventAuxInfo& eventAuxInfo )
252			{
253			std::stringstream eventStringStream;
254
255			std::vector<VHbbEventAuxInfo::ParticleMCInfo> zMCInfo=eventAuxInfo.mcZ;
256			std::vector<VHbbEventAuxInfo::ParticleMCInfo> higgsMCInfo=eventAuxInfo.mcH;
257
258			for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator higgsInfo=higgsMCInfo.begin(); higgsInfo!=higgsMCInfo.end(); ++higgsInfo )
259			{
260			if( higgsInfo->dauid.size()==1 && std::abs(higgsInfo->dauid[0])==25 ) continue;
261
262			eventStringStream << "H->(";
263			for( std::vector<int>::const_iterator daughterID=higgsInfo->dauid.begin(); daughterID!=higgsInfo->dauid.end(); ++daughterID )
264			{
265			if( daughterID!=higgsInfo->dauid.begin() ) eventStringStream << ",";
266			eventStringStream << *daughterID;
267			}
268			eventStringStream << ") ";
269			}
270
271			for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator zInfo=zMCInfo.begin(); zInfo!=zMCInfo.end(); ++zInfo )
272			{
273			if( zInfo->dauid.empty() ) continue;
274			if( zInfo->dauid.size()==1 && std::abs(zInfo->dauid[0])==23 ) continue;
275			if( zInfo->dauid.size()==2 && std::abs(zInfo->dauid[0])==23 && std::abs(zInfo->dauid[1])==25 ) continue;
276
277			eventStringStream << "Z->(";
278			for( std::vector<int>::const_iterator daughterID=zInfo->dauid.begin(); daughterID!=zInfo->dauid.end(); ++daughterID )
279			{
280			if( *daughterID==23 ) continue;
281			if( zInfo->dauid.size()>2 && *daughterID==22 ) continue;
282			if( daughterID!=zInfo->dauid.begin() ) eventStringStream << ",";
283			eventStringStream << *daughterID;
284			}
285			eventStringStream << ") ";
286			}
287
288	grimes	1.2	if( eventStringStream.str().empty() ) return "otherBackground";
289			else return eventStringStream.str();
290	grimes	1.1	}
291