ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/VHbbAnalysisCode/src/LightweightAnalyser.cpp
Revision: 1.4
Committed: Fri Apr 27 13:52:14 2012 UTC (13 years ago) by grimes
Branch: MAIN
Changes since 1.3: +8 -4 lines
Log Message:
Added code to analyse muon isolation.

File Contents

# 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.1
8     #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
9     #include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc" // Not entirely sure why this is included and not linked to
10     #include "DataFormats/FWLite/interface/Event.h"
11     #include "DataFormats/FWLite/interface/Handle.h"
12     #include "DataFormats/PatCandidates/interface/Muon.h"
13    
14    
15     #include <TFile.h>
16     #include <TH1.h>
17     #include <TH2.h>
18    
19    
20     const bool trkupgradeanalysis::LightweightAnalyser::verbose_=false;
21     const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdZ_=20;
22     const double trkupgradeanalysis::LightweightAnalyser::jetPtThresholdW_=20;
23     const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsZ_=false;
24     const bool trkupgradeanalysis::LightweightAnalyser::useHighestPtHiggsW_=true;
25    
26     namespace // Use the unnamed namespace
27     {
28     TDirectory* createDirectory( const std::string& fullPath, TDirectory* pParent )
29     {
30     if( pParent==NULL ) throw std::runtime_error( "The parent directory is a Null pointer" );
31    
32     TDirectory* pSubDirectory=pParent;
33     size_t currentPosition=0;
34     size_t nextSlash;
35     do
36     {
37     nextSlash=fullPath.find_first_of('/', currentPosition );
38     std::string directoryName=fullPath.substr(currentPosition,nextSlash-currentPosition);
39     currentPosition=nextSlash+1;
40    
41     TDirectory* pNextSubDirectory=pSubDirectory->GetDirectory( directoryName.c_str() );
42     if( pNextSubDirectory==NULL ) pNextSubDirectory=pSubDirectory->mkdir( directoryName.c_str() );
43     if( pNextSubDirectory==NULL ) throw std::runtime_error( "Couldn't create the root directory \""+directoryName+"\"" );
44     pSubDirectory=pNextSubDirectory;
45    
46     } while( nextSlash!=std::string::npos );
47    
48     return pSubDirectory;
49     }
50     }
51    
52     trkupgradeanalysis::LightweightAnalyser::LightweightAnalyser( const std::string& outputFilename )
53     : zFinder_( verbose_, jetPtThresholdZ_, useHighestPtHiggsZ_ ),
54     wFinder_( verbose_, jetPtThresholdW_, useHighestPtHiggsW_ ),
55     numberOfEvents_(0),
56     numberOfEventsWithAtLeastOneZOrWCandidate_(0),
57     pPlotsForEverything_(NULL)
58     {
59     // All of these are empty smart pointers, so I need to call reset on them
60     pOutputFile_.reset( new TFile( outputFilename.c_str(), "RECREATE" ) );
61    
62 grimes 1.2
63     using namespace trkupgradeanalysis; // Don't really need this but I like to be explicit (ooh err..)
64     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZee(120)) );
65     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(120)) );
66     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionWen(120)) );
67     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionWmun(120)) );
68     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumuWithoutAdditionalJetsCut(120)) );
69     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(110)) );
70     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(115)) );
71     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(125)) );
72     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(130)) );
73     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(135)) );
74 grimes 1.3 cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumu(110-20,135+10)) );
75     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new SignalSelectionZmumuWithCSVCutsSwitched(110-20,135+10)) );
76 grimes 1.2
77     //
78     // These are only temporary while I'm studying the number of additional jets
79     //
80     using namespace trkupgradeanalysis::additionaljetstudy;
81     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide2() ) );
82     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide3Plot1() ) );
83     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide3Plot2() ) );
84     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide3Plot3() ) );
85     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide4Plot2() ) );
86     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide5Plot1() ) );
87     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide5Plot2() ) );
88     cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>( new WilkenSlide5Plot2AssumingTypo() ) );
89    
90     //
91     // Don't need the following because I'm not studying the background estimation regions
92     //
93    
94     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VlightRegionHWmun) );
95     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VlightRegionHWen) );
96     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VlightRegionHZmumu) );
97     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VlightRegionHZee) );
98     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new TTbarRegionHWmun) );
99     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new TTbarRegionHWen) );
100     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new TTbarRegionHZmumu) );
101     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new TTbarRegionHZee) );
102     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VbbRegionHWmun) );
103     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VbbRegionHWen) );
104     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VbbRegionHZmumu) );
105     //cutsToApply_.push_back( boost::shared_ptr<VHbbCandidateCutSet>(new VbbRegionHZee) );
106 grimes 1.1 }
107    
108     trkupgradeanalysis::LightweightAnalyser::~LightweightAnalyser()
109     {
110     save();
111     pOutputFile_->Close();
112     }
113    
114     void trkupgradeanalysis::LightweightAnalyser::changeDirectory( const std::string& directory )
115     {
116     std::map<std::string,PlotsSplitByMCType>::iterator iPlots=plotsForEventTypes_.find(directory);
117    
118     if( iPlots==plotsForEventTypes_.end() )
119     {
120     pPlotsForEverything_=&plotsForEventTypes_[directory]; // This will create a new entry in the map
121    
122     // Only book the histograms for all MC event types now. I don't know what MC event types there
123     // until I look at the events so I need to wait until then
124     bookHistograms( pPlotsForEverything_->allEventTypes_, directory+"/allEvents" );
125     }
126     else pPlotsForEverything_=&iPlots->second; // This directory has already been created and the histograms booked, so reuse those
127    
128     currentDirectory_=directory;
129     }
130    
131     void trkupgradeanalysis::LightweightAnalyser::save()
132     {
133     pOutputFile_->Write(0,TObject::kOverwrite);
134     }
135    
136     void trkupgradeanalysis::LightweightAnalyser::bookHistograms( trkupgradeanalysis::LightweightAnalyser::AllPlots& plots, const std::string& directoryName )
137     {
138     TDirectory* pRootDirectory=pOutputFile_->GetDirectory("/");
139     TDirectory* pDirectory=createDirectory( directoryName, pRootDirectory );
140     // TDirectory* pDirectory=pRootDirectory->GetDirectory( directoryName.c_str() );
141     // if( pDirectory==NULL ) pDirectory=pRootDirectory->mkdir( directoryName.c_str() );
142    
143     TDirectory* pSubDirectory=pDirectory->mkdir( "VHbbCandidates" );
144    
145     plots.pNumberOfCandidates_=new TH2F( "numberOfCandidates","Number of candidates",5,-0.5,4.5,5,-0.5,4.5);
146     plots.pNumberOfCandidates_->SetDirectory(pSubDirectory);
147     plots.candidateHistogramsForAllEvents_.book( pSubDirectory );
148     plots.pNumberOfCandidates_->GetXaxis()->SetTitle( "Z candidates" );
149     plots.pNumberOfCandidates_->GetYaxis()->SetTitle( "W candidates" );
150    
151     pSubDirectory=pDirectory->mkdir( "MuonInfos" );
152     plots.allMuons_.book( pSubDirectory );
153    
154     pSubDirectory=pDirectory->mkdir( "ElectronInfos" );
155     plots.allElectrons_.book( pSubDirectory );
156    
157 grimes 1.4 pSubDirectory=pDirectory->mkdir( "IsolationStudy" );
158     plots.isolationStudy_.book( pSubDirectory );
159    
160 grimes 1.1 pSubDirectory=pDirectory->mkdir( "MonteCarloTruth" );
161     plots.monteCarloInfo_.book( pSubDirectory );
162    
163     pSubDirectory=pDirectory->mkdir( "Cuts" );
164    
165     for( std::vector< boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet> >::const_iterator iCutSet=cutsToApply_.begin();
166     iCutSet!=cutsToApply_.end(); ++iCutSet )
167     {
168     TDirectory* pSubSubDirectory=pSubDirectory->mkdir( (*iCutSet)->name().c_str() );
169    
170     CutSetPlotSet newPlotSet( *iCutSet );
171     newPlotSet.book( pSubSubDirectory );
172     plots.cutCollectionPlotSets_.push_back( newPlotSet );
173     }
174     }
175    
176 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 )
177 grimes 1.1 {
178     pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo );
179     pPlotsToFill->allElectrons_.fill( vhbbEvent.eleInfo );
180 grimes 1.4 if( pVHbbEventAuxInfo ) pPlotsToFill->monteCarloInfo_.fill( *pVHbbEventAuxInfo );
181     pPlotsToFill->isolationStudy_.fill( vhbbEvent, pVHbbEventAuxInfo ); // This class checks if is null itself
182 grimes 1.1
183     for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
184     {
185     pPlotsToFill->candidateHistogramsForAllEvents_.fill( *iCandidate );
186    
187     // Loop over all of the cuts. If this candidate passes the cuts then fill the plots
188     for( std::vector<CutSetPlotSet>::iterator iCutPlotSet=pPlotsToFill->cutCollectionPlotSets_.begin();
189     iCutPlotSet!=pPlotsToFill->cutCollectionPlotSets_.end(); ++iCutPlotSet )
190     {
191     iCutPlotSet->fill( *iCandidate );
192     }
193     }
194    
195     pPlotsToFill->pNumberOfCandidates_->Fill( zCandidates.size(), wCandidates.size() );
196     }
197    
198     void trkupgradeanalysis::LightweightAnalyser::processEvent( const fwlite::Event& event )
199     {
200     std::vector<VHbbCandidate> zCandidates;
201     std::vector<VHbbCandidate> wCandidates;
202    
203     fwlite::Handle<VHbbEvent> vhbbHandle;
204     vhbbHandle.getByLabel( event, "HbbAnalyzerNew" );
205     const VHbbEvent* pVHbbEvent=vhbbHandle.product();
206     if( pVHbbEvent == 0 )
207     {
208     std::cerr << "Couldn't get the VHbbEvent" << std::endl;
209     return;
210     }
211    
212     AllPlots* pCurrentPlotsForThisEventType=NULL;
213     fwlite::Handle<VHbbEventAuxInfo> vhbbAuxInfoHandle;
214     vhbbAuxInfoHandle.getByLabel( event, "HbbAnalyzerNew" );
215     const VHbbEventAuxInfo* pVHbbEventAuxInfo=vhbbAuxInfoHandle.product();
216 grimes 1.2 if( pVHbbEventAuxInfo == 0 ) std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
217 grimes 1.1 else
218     {
219     std::string eventType=eventTypeFromMC( *pVHbbEventAuxInfo );
220    
221     std::map<std::string,AllPlots>::iterator iPlotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_.find(eventType);
222    
223     if( iPlotsForThisEventType==pPlotsForEverything_->mcEventTypePlots_.end() )
224     {
225     // This event type hasn't been encountered before so I need to book the histograms
226     AllPlots& plotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_[eventType];
227     bookHistograms( plotsForThisEventType, currentDirectory_+"/"+eventType );
228     pCurrentPlotsForThisEventType=&plotsForThisEventType;
229     }
230     else pCurrentPlotsForThisEventType=&(iPlotsForThisEventType->second);
231     }
232    
233     ++numberOfEvents_;
234    
235     zFinder_.run( pVHbbEvent, zCandidates );
236     wFinder_.run( pVHbbEvent, wCandidates );
237    
238 grimes 1.4 fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
239     if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, pVHbbEventAuxInfo, zCandidates, wCandidates );
240 grimes 1.1
241    
242     if( !zCandidates.empty() || !wCandidates.empty() ) ++numberOfEventsWithAtLeastOneZOrWCandidate_;
243    
244     }
245    
246     void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile )
247     {
248     if( pPlotsForEverything_==NULL ) changeDirectory("allEvents");
249    
250     fwlite::Event event( pInputFile );
251     for( event.toBegin(); !event.atEnd(); ++event )
252     {
253     processEvent( event );
254     }
255    
256     std::cout << "There were " << numberOfEventsWithAtLeastOneZOrWCandidate_ << " events with a Z or W candidate out of " << numberOfEvents_ << " events." << std::endl;
257     }
258    
259     std::string trkupgradeanalysis::LightweightAnalyser::eventTypeFromMC( const VHbbEventAuxInfo& eventAuxInfo )
260     {
261     std::stringstream eventStringStream;
262    
263     std::vector<VHbbEventAuxInfo::ParticleMCInfo> zMCInfo=eventAuxInfo.mcZ;
264     std::vector<VHbbEventAuxInfo::ParticleMCInfo> higgsMCInfo=eventAuxInfo.mcH;
265    
266     for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator higgsInfo=higgsMCInfo.begin(); higgsInfo!=higgsMCInfo.end(); ++higgsInfo )
267     {
268     if( higgsInfo->dauid.size()==1 && std::abs(higgsInfo->dauid[0])==25 ) continue;
269    
270     eventStringStream << "H->(";
271     for( std::vector<int>::const_iterator daughterID=higgsInfo->dauid.begin(); daughterID!=higgsInfo->dauid.end(); ++daughterID )
272     {
273     if( daughterID!=higgsInfo->dauid.begin() ) eventStringStream << ",";
274     eventStringStream << *daughterID;
275     }
276     eventStringStream << ") ";
277     }
278    
279     for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator zInfo=zMCInfo.begin(); zInfo!=zMCInfo.end(); ++zInfo )
280     {
281     if( zInfo->dauid.empty() ) continue;
282     if( zInfo->dauid.size()==1 && std::abs(zInfo->dauid[0])==23 ) continue;
283     if( zInfo->dauid.size()==2 && std::abs(zInfo->dauid[0])==23 && std::abs(zInfo->dauid[1])==25 ) continue;
284    
285     eventStringStream << "Z->(";
286     for( std::vector<int>::const_iterator daughterID=zInfo->dauid.begin(); daughterID!=zInfo->dauid.end(); ++daughterID )
287     {
288     if( *daughterID==23 ) continue;
289     if( zInfo->dauid.size()>2 && *daughterID==22 ) continue;
290     if( daughterID!=zInfo->dauid.begin() ) eventStringStream << ",";
291     eventStringStream << *daughterID;
292     }
293     eventStringStream << ") ";
294     }
295    
296 grimes 1.2 if( eventStringStream.str().empty() ) return "otherBackground";
297     else return eventStringStream.str();
298 grimes 1.1 }
299