| 1 |
#include "TrkUpgradeAnalysis/VHbb/interface/MCInfoPlotSet.h"
|
| 2 |
|
| 3 |
#include <stdexcept>
|
| 4 |
#include <sstream>
|
| 5 |
|
| 6 |
#include <TDirectory.h>
|
| 7 |
#include <TH1F.h>
|
| 8 |
|
| 9 |
|
| 10 |
trkupgradeanalysis::MCInfoPlotSet::MCInfoPlotSet()
|
| 11 |
: histogramHaveBeenBooked_(false)
|
| 12 |
{
|
| 13 |
// No operation besides the initialiser list.
|
| 14 |
}
|
| 15 |
|
| 16 |
void trkupgradeanalysis::MCInfoPlotSet::book( TDirectory* pDirectory )
|
| 17 |
{
|
| 18 |
if( histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MCInfoPlotSet::book() - histograms have already been booked" );
|
| 19 |
|
| 20 |
//
|
| 21 |
// Note that the root file which TDirectory is part of takes ownership of all
|
| 22 |
// of these objects, so I don't need to (and shouldn't) delete them when I'm
|
| 23 |
// finished.
|
| 24 |
//
|
| 25 |
|
| 26 |
|
| 27 |
pEventString_=new TH1F( "mcEventType", "Event type from MC", 10, 0.5, 10.5 );
|
| 28 |
pEventString_->SetDirectory(pDirectory);
|
| 29 |
|
| 30 |
pNumberOfPrimaryVertices_=new TH1F( "numberOfPrimaryVertices", "Number of primary vertices", 71, -0.5, 70.5 );
|
| 31 |
pNumberOfPrimaryVertices_->SetDirectory(pDirectory);
|
| 32 |
|
| 33 |
pNumberOfBunchCrossings_=new TH1F( "numberOfBunchCrossings", "Number of bunch crossings", 71, -0.5, 70.5 );
|
| 34 |
pNumberOfInteractionsPerBunchCrossing_=new TH1F( "numberOfInteractionsPerBunchCrossing", "Number of interactions per bunch crossing", 71, -0.5, 70.5 );
|
| 35 |
pTotalInteractionsPerEvent_=new TH1F( "totalInteractionsPerEvent", "Total number of interactions per event", 71, -0.5, 70.5 );
|
| 36 |
pNumberOfBunchCrossings_->SetDirectory(pDirectory);
|
| 37 |
pNumberOfInteractionsPerBunchCrossing_->SetDirectory(pDirectory);
|
| 38 |
pTotalInteractionsPerEvent_->SetDirectory(pDirectory);
|
| 39 |
|
| 40 |
histogramHaveBeenBooked_=true;
|
| 41 |
}
|
| 42 |
|
| 43 |
void trkupgradeanalysis::MCInfoPlotSet::fill( const VHbbEventAuxInfo& eventAuxInfo )
|
| 44 |
{
|
| 45 |
if( !histogramHaveBeenBooked_ ) throw std::runtime_error( "trkupgradeanalysis::MCInfoPlotSet::book() - histograms have not been booked" );
|
| 46 |
|
| 47 |
std::stringstream eventStringStream;
|
| 48 |
|
| 49 |
std::vector<VHbbEventAuxInfo::ParticleMCInfo> zMCInfo=eventAuxInfo.mcZ;
|
| 50 |
std::vector<VHbbEventAuxInfo::ParticleMCInfo> higgsMCInfo=eventAuxInfo.mcH;
|
| 51 |
|
| 52 |
for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator higgsInfo=higgsMCInfo.begin(); higgsInfo!=higgsMCInfo.end(); ++higgsInfo )
|
| 53 |
{
|
| 54 |
if( higgsInfo->dauid.size()==1 && std::abs(higgsInfo->dauid[0])==25 ) continue;
|
| 55 |
|
| 56 |
eventStringStream << "H->(";
|
| 57 |
for( std::vector<int>::const_iterator daughterID=higgsInfo->dauid.begin(); daughterID!=higgsInfo->dauid.end(); ++daughterID )
|
| 58 |
{
|
| 59 |
if( daughterID!=higgsInfo->dauid.begin() ) eventStringStream << ",";
|
| 60 |
eventStringStream << *daughterID;
|
| 61 |
}
|
| 62 |
eventStringStream << ") ";
|
| 63 |
}
|
| 64 |
|
| 65 |
for( std::vector<VHbbEventAuxInfo::ParticleMCInfo>::const_iterator zInfo=zMCInfo.begin(); zInfo!=zMCInfo.end(); ++zInfo )
|
| 66 |
{
|
| 67 |
if( zInfo->dauid.empty() ) continue;
|
| 68 |
if( zInfo->dauid.size()==1 && std::abs(zInfo->dauid[0])==23 ) continue;
|
| 69 |
if( zInfo->dauid.size()==2 && std::abs(zInfo->dauid[0])==23 && std::abs(zInfo->dauid[1])==25 ) continue;
|
| 70 |
|
| 71 |
eventStringStream << "Z->(";
|
| 72 |
for( std::vector<int>::const_iterator daughterID=zInfo->dauid.begin(); daughterID!=zInfo->dauid.end(); ++daughterID )
|
| 73 |
{
|
| 74 |
if( *daughterID==23 ) continue;
|
| 75 |
if( zInfo->dauid.size()>2 && *daughterID==22 ) continue;
|
| 76 |
if( daughterID!=zInfo->dauid.begin() ) eventStringStream << ",";
|
| 77 |
eventStringStream << *daughterID;
|
| 78 |
}
|
| 79 |
eventStringStream << ") ";
|
| 80 |
}
|
| 81 |
|
| 82 |
// Now that I've created the string describing the MC truth type of the event, see
|
| 83 |
// if it's already been used. If so find out which bin it's being plotted in. Note that
|
| 84 |
// "binNumber" here isn't strictly speaking the bin number, it's the x-coord used to fill
|
| 85 |
// the correct bin.
|
| 86 |
int binNumber=0;
|
| 87 |
std::map<std::string,int>::const_iterator iBinNumber=stringToBinMap_.find( eventStringStream.str() );
|
| 88 |
if( iBinNumber==stringToBinMap_.end() )
|
| 89 |
{
|
| 90 |
// hasn't been plotted before, so I need to create an entry in the map
|
| 91 |
// and use the next available bin.
|
| 92 |
binNumber=stringToBinMap_.size()+1;
|
| 93 |
stringToBinMap_[eventStringStream.str()]=binNumber;
|
| 94 |
|
| 95 |
// Give the histogram bin the correct label
|
| 96 |
pEventString_->GetXaxis()->SetBinLabel( binNumber, eventStringStream.str().c_str() );
|
| 97 |
}
|
| 98 |
else binNumber=iBinNumber->second;
|
| 99 |
|
| 100 |
pEventString_->Fill( binNumber );
|
| 101 |
|
| 102 |
pNumberOfPrimaryVertices_->Fill( eventAuxInfo.pvInfo.nVertices );
|
| 103 |
|
| 104 |
// Loop over the pile up data
|
| 105 |
unsigned int totalNumberOfInteractions=0;
|
| 106 |
unsigned int numberOfBunchCrossings=0;
|
| 107 |
for( std::map<int,unsigned int>::const_iterator iBXInteractionPair=eventAuxInfo.puInfo.pus.begin(); iBXInteractionPair!=eventAuxInfo.puInfo.pus.end(); ++iBXInteractionPair )
|
| 108 |
{
|
| 109 |
const unsigned int& numberOfInteractions=iBXInteractionPair->second;
|
| 110 |
totalNumberOfInteractions+=numberOfInteractions;
|
| 111 |
pNumberOfInteractionsPerBunchCrossing_->Fill( numberOfInteractions );
|
| 112 |
++numberOfBunchCrossings;
|
| 113 |
}
|
| 114 |
pNumberOfBunchCrossings_->Fill( numberOfBunchCrossings );
|
| 115 |
pTotalInteractionsPerEvent_->Fill( totalNumberOfInteractions );
|
| 116 |
|
| 117 |
}
|
| 118 |
|
