ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/VHbbAnalysisCode/src/LightweightAnalyser.cpp
Revision: 1.3
Committed: Tue Mar 27 23:50:44 2012 UTC (13 years, 1 month ago) by grimes
Branch: MAIN
Changes since 1.2: +2 -0 lines
Log Message:
Various improvements.

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     pSubDirectory=pDirectory->mkdir( "MonteCarloTruth" );
158     plots.monteCarloInfo_.book( pSubDirectory );
159    
160     pSubDirectory=pDirectory->mkdir( "Cuts" );
161    
162     for( std::vector< boost::shared_ptr<trkupgradeanalysis::VHbbCandidateCutSet> >::const_iterator iCutSet=cutsToApply_.begin();
163     iCutSet!=cutsToApply_.end(); ++iCutSet )
164     {
165     TDirectory* pSubSubDirectory=pSubDirectory->mkdir( (*iCutSet)->name().c_str() );
166    
167     CutSetPlotSet newPlotSet( *iCutSet );
168     newPlotSet.book( pSubSubDirectory );
169     plots.cutCollectionPlotSets_.push_back( newPlotSet );
170     }
171     }
172    
173     void trkupgradeanalysis::LightweightAnalyser::fillAllPlotsStructure( AllPlots* pPlotsToFill, const VHbbEvent& vhbbEvent, const std::vector<VHbbCandidate>& zCandidates, const std::vector<VHbbCandidate>& wCandidates )
174     {
175     pPlotsToFill->allMuons_.fill( vhbbEvent.muInfo );
176     pPlotsToFill->allElectrons_.fill( vhbbEvent.eleInfo );
177    
178     for( std::vector<VHbbCandidate>::const_iterator iCandidate=zCandidates.begin(); iCandidate!=zCandidates.end(); ++iCandidate )
179     {
180     pPlotsToFill->candidateHistogramsForAllEvents_.fill( *iCandidate );
181    
182     // Loop over all of the cuts. If this candidate passes the cuts then fill the plots
183     for( std::vector<CutSetPlotSet>::iterator iCutPlotSet=pPlotsToFill->cutCollectionPlotSets_.begin();
184     iCutPlotSet!=pPlotsToFill->cutCollectionPlotSets_.end(); ++iCutPlotSet )
185     {
186     iCutPlotSet->fill( *iCandidate );
187     }
188     }
189    
190     pPlotsToFill->pNumberOfCandidates_->Fill( zCandidates.size(), wCandidates.size() );
191     }
192    
193     void trkupgradeanalysis::LightweightAnalyser::processEvent( const fwlite::Event& event )
194     {
195     std::vector<VHbbCandidate> zCandidates;
196     std::vector<VHbbCandidate> wCandidates;
197    
198     fwlite::Handle<VHbbEvent> vhbbHandle;
199     vhbbHandle.getByLabel( event, "HbbAnalyzerNew" );
200     const VHbbEvent* pVHbbEvent=vhbbHandle.product();
201     if( pVHbbEvent == 0 )
202     {
203     std::cerr << "Couldn't get the VHbbEvent" << std::endl;
204     return;
205     }
206    
207     AllPlots* pCurrentPlotsForThisEventType=NULL;
208     fwlite::Handle<VHbbEventAuxInfo> vhbbAuxInfoHandle;
209     vhbbAuxInfoHandle.getByLabel( event, "HbbAnalyzerNew" );
210     const VHbbEventAuxInfo* pVHbbEventAuxInfo=vhbbAuxInfoHandle.product();
211 grimes 1.2 if( pVHbbEventAuxInfo == 0 ) std::cerr << "Couldn't get the VHbbEventAuxInfo" << std::endl;
212 grimes 1.1 else
213     {
214     pPlotsForEverything_->allEventTypes_.monteCarloInfo_.fill( *pVHbbEventAuxInfo );
215     std::string eventType=eventTypeFromMC( *pVHbbEventAuxInfo );
216    
217     std::map<std::string,AllPlots>::iterator iPlotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_.find(eventType);
218    
219     if( iPlotsForThisEventType==pPlotsForEverything_->mcEventTypePlots_.end() )
220     {
221     // This event type hasn't been encountered before so I need to book the histograms
222     AllPlots& plotsForThisEventType=pPlotsForEverything_->mcEventTypePlots_[eventType];
223     bookHistograms( plotsForThisEventType, currentDirectory_+"/"+eventType );
224     pCurrentPlotsForThisEventType=&plotsForThisEventType;
225     }
226     else pCurrentPlotsForThisEventType=&(iPlotsForThisEventType->second);
227     }
228    
229     ++numberOfEvents_;
230    
231     zFinder_.run( pVHbbEvent, zCandidates );
232     wFinder_.run( pVHbbEvent, wCandidates );
233    
234     fillAllPlotsStructure( &pPlotsForEverything_->allEventTypes_, *pVHbbEvent, zCandidates, wCandidates );
235     if( pCurrentPlotsForThisEventType!=NULL ) fillAllPlotsStructure( pCurrentPlotsForThisEventType, *pVHbbEvent, zCandidates, wCandidates );
236    
237    
238     if( !zCandidates.empty() || !wCandidates.empty() ) ++numberOfEventsWithAtLeastOneZOrWCandidate_;
239    
240     }
241    
242     void trkupgradeanalysis::LightweightAnalyser::processFile( TFile* pInputFile )
243     {
244     if( pPlotsForEverything_==NULL ) changeDirectory("allEvents");
245    
246     fwlite::Event event( pInputFile );
247     for( event.toBegin(); !event.atEnd(); ++event )
248     {
249     processEvent( event );
250     }
251    
252     std::cout << "There were " << numberOfEventsWithAtLeastOneZOrWCandidate_ << " events with a Z or W candidate out of " << numberOfEvents_ << " events." << std::endl;
253     }
254    
255     std::string trkupgradeanalysis::LightweightAnalyser::eventTypeFromMC( const VHbbEventAuxInfo& eventAuxInfo )
256     {
257     std::stringstream eventStringStream;
258    
259     std::vector<VHbbEventAuxInfo::ParticleMCInfo> zMCInfo=eventAuxInfo.mcZ;
260     std::vector<VHbbEventAuxInfo::ParticleMCInfo> higgsMCInfo=eventAuxInfo.mcH;
261    
262     for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator higgsInfo=higgsMCInfo.begin(); higgsInfo!=higgsMCInfo.end(); ++higgsInfo )
263     {
264     if( higgsInfo->dauid.size()==1 && std::abs(higgsInfo->dauid[0])==25 ) continue;
265    
266     eventStringStream << "H->(";
267     for( std::vector<int>::const_iterator daughterID=higgsInfo->dauid.begin(); daughterID!=higgsInfo->dauid.end(); ++daughterID )
268     {
269     if( daughterID!=higgsInfo->dauid.begin() ) eventStringStream << ",";
270     eventStringStream << *daughterID;
271     }
272     eventStringStream << ") ";
273     }
274    
275     for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator zInfo=zMCInfo.begin(); zInfo!=zMCInfo.end(); ++zInfo )
276     {
277     if( zInfo->dauid.empty() ) continue;
278     if( zInfo->dauid.size()==1 && std::abs(zInfo->dauid[0])==23 ) continue;
279     if( zInfo->dauid.size()==2 && std::abs(zInfo->dauid[0])==23 && std::abs(zInfo->dauid[1])==25 ) continue;
280    
281     eventStringStream << "Z->(";
282     for( std::vector<int>::const_iterator daughterID=zInfo->dauid.begin(); daughterID!=zInfo->dauid.end(); ++daughterID )
283     {
284     if( *daughterID==23 ) continue;
285     if( zInfo->dauid.size()>2 && *daughterID==22 ) continue;
286     if( daughterID!=zInfo->dauid.begin() ) eventStringStream << ",";
287     eventStringStream << *daughterID;
288     }
289     eventStringStream << ") ";
290     }
291    
292 grimes 1.2 if( eventStringStream.str().empty() ) return "otherBackground";
293     else return eventStringStream.str();
294 grimes 1.1 }
295