| 1 |
#include "OSUT3Analysis/AnaTools/interface/OSUAnalysis.h"
|
| 2 |
|
| 3 |
OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
|
| 4 |
// Retrieve parameters from the configuration file.
|
| 5 |
muons_ (cfg.getParameter<edm::InputTag> ("muons")),
|
| 6 |
electrons_ (cfg.getParameter<edm::InputTag> ("electrons")),
|
| 7 |
channels_ (cfg.getParameter<vector<edm::ParameterSet> >("channels"))
|
| 8 |
|
| 9 |
{
|
| 10 |
TH1::SetDefaultSumw2 ();
|
| 11 |
|
| 12 |
// Construct Cutflow Objects. These store the results of cut decisions and
|
| 13 |
// handle filling cut flow histograms.
|
| 14 |
masterCutFlow_ = new CutFlow (fs_);
|
| 15 |
std::vector<TFileDirectory> directories;
|
| 16 |
|
| 17 |
|
| 18 |
//loop over all channels (event selections)
|
| 19 |
for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
|
| 20 |
|
| 21 |
//get name of channel
|
| 22 |
string channelName (channels_.at(currentChannel).getParameter<string>("name"));
|
| 23 |
channels.push_back(channelName);
|
| 24 |
TString channelLabel = channelName;
|
| 25 |
|
| 26 |
|
| 27 |
//create cutFlow for this channel
|
| 28 |
cutFlows_.push_back (new CutFlow (fs_, channelName));
|
| 29 |
|
| 30 |
//book a directory in the output file with the name of the channel
|
| 31 |
TFileDirectory subDir = fs_->mkdir( channelName );
|
| 32 |
directories.push_back(subDir);
|
| 33 |
std::map<std::string, TH1D*> histoMap;
|
| 34 |
oneDHists_.push_back(histoMap);
|
| 35 |
|
| 36 |
|
| 37 |
//muon histograms
|
| 38 |
oneDHists_.at(currentChannel)["muonPt"] = directories.at(currentChannel).make<TH1D> ("muonPt",channelLabel+" channel: Muon Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
|
| 39 |
oneDHists_.at(currentChannel)["muonEta"] = directories.at(currentChannel).make<TH1D> ("muonEta",channelLabel+" channel: Muon Eta; #eta", 100, -5, 5);
|
| 40 |
oneDHists_.at(currentChannel)["muonD0"] = directories.at(currentChannel).make<TH1D> ("muonD0",channelLabel+" channel: Muon d_{0}; d_{0} [cm] ", 1000, -1, 1);
|
| 41 |
oneDHists_.at(currentChannel)["muonAbsD0"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0",channelLabel+" channel: Muon d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
|
| 42 |
oneDHists_.at(currentChannel)["muonD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonD0Sig",channelLabel+" channel: Muon d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
|
| 43 |
oneDHists_.at(currentChannel)["muonAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0Sig",channelLabel+" channel: Muon d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
|
| 44 |
oneDHists_.at(currentChannel)["muonIso"] = directories.at(currentChannel).make<TH1D> ("muonIso",channelLabel+" channel: Muon Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
|
| 45 |
oneDHists_.at(currentChannel)["numMuons"] = directories.at(currentChannel).make<TH1D> ("numMuons",channelLabel+" channel: Number of Selected Muons; # muons", 10, 0, 10);
|
| 46 |
|
| 47 |
|
| 48 |
//electron histograms
|
| 49 |
oneDHists_.at(currentChannel)["electronPt"] = directories.at(currentChannel).make<TH1D> ("electronPt",channelLabel+" channel: Electron Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
|
| 50 |
oneDHists_.at(currentChannel)["electronEta"] = directories.at(currentChannel).make<TH1D> ("electronEta",channelLabel+" channel: Electron Eta; #eta", 50, -5, 5);
|
| 51 |
oneDHists_.at(currentChannel)["electronD0"] = directories.at(currentChannel).make<TH1D> ("electronD0",channelLabel+" channel: Electron d_{0}; d_{0} [cm] ", 1000, -1, 1);
|
| 52 |
oneDHists_.at(currentChannel)["electronAbsD0"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0",channelLabel+" channel: Electron d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
|
| 53 |
oneDHists_.at(currentChannel)["electronD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronD0Sig",channelLabel+" channel: Electron d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
|
| 54 |
oneDHists_.at(currentChannel)["electronAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0Sig",channelLabel+" channel: Electron d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
|
| 55 |
oneDHists_.at(currentChannel)["electronIso"] = directories.at(currentChannel).make<TH1D> ("electronIso",channelLabel+" channel: Electron Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
|
| 56 |
oneDHists_.at(currentChannel)["electronFbrem"] = directories.at(currentChannel).make<TH1D> ("electronFbrem",channelLabel+" channel: Electron Brem. Energy Fraction; fbrem", 1000, 0, 2);
|
| 57 |
oneDHists_.at(currentChannel)["numElectrons"] = directories.at(currentChannel).make<TH1D> ("numElectrons",channelLabel+" channel: Number of Selected Electrons; # electrons", 10, 0, 10);
|
| 58 |
|
| 59 |
|
| 60 |
|
| 61 |
//get list of cuts for this channel
|
| 62 |
vector<edm::ParameterSet> cuts_ (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
|
| 63 |
|
| 64 |
//loop over and parse all cuts
|
| 65 |
for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
|
| 66 |
|
| 67 |
//store input collection for cut
|
| 68 |
string inputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
|
| 69 |
inputCollections[channelName].push_back(inputCollection);
|
| 70 |
|
| 71 |
//split cut string into parts and store them
|
| 72 |
string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
|
| 73 |
std::vector<string> cutStringVector = splitString(cutString);
|
| 74 |
cutVariables[channelName].push_back(cutStringVector.at(0)); // variable to cut on
|
| 75 |
cutComparativeOperators[channelName].push_back(cutStringVector.at(1)); // comparison to make
|
| 76 |
cutValues[channelName].push_back(atof(cutStringVector.at(2).c_str())); // threshold value to pass cut
|
| 77 |
|
| 78 |
//get number of objects required to pass cut for event to pass
|
| 79 |
string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
|
| 80 |
std::vector<string> numberRequiredVector = splitString(numberRequiredString);
|
| 81 |
|
| 82 |
// determine number required if comparison contains "="
|
| 83 |
int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
|
| 84 |
if(numberRequiredVector.at(0) == ">") numberRequiredInt++;
|
| 85 |
else if(numberRequiredVector.at(0) == "<") numberRequiredInt--;
|
| 86 |
|
| 87 |
|
| 88 |
numberRequiredInts[channelName].push_back(numberRequiredInt); // number of objects required to pass the cut
|
| 89 |
eventComparativeOperators[channelName].push_back(numberRequiredVector.at(0)); // comparison to make
|
| 90 |
|
| 91 |
|
| 92 |
if(cuts_.at(currentCut).exists("alias")){
|
| 93 |
cutNames[channelName].push_back(cuts_.at(currentCut).getParameter<string> ("alias"));
|
| 94 |
}
|
| 95 |
else{
|
| 96 |
//construct string for cutflow table
|
| 97 |
bool plural = numberRequiredInt != 1;
|
| 98 |
string collectionString = plural ? inputCollection : inputCollection.substr(0, inputCollection.size()-1);
|
| 99 |
string cutName = numberRequiredString + " " + collectionString + " with " + cutString;
|
| 100 |
cutNames[channelName].push_back(cutName);
|
| 101 |
}
|
| 102 |
|
| 103 |
}//end loop over cuts
|
| 104 |
|
| 105 |
|
| 106 |
|
| 107 |
}//end loop over channels
|
| 108 |
|
| 109 |
|
| 110 |
|
| 111 |
|
| 112 |
}
|
| 113 |
|
| 114 |
OSUAnalysis::~OSUAnalysis ()
|
| 115 |
{
|
| 116 |
// Destroying the CutFlow objects causes the cut flow numbers and time
|
| 117 |
// information to be printed to standard output.
|
| 118 |
for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
|
| 119 |
delete cutFlows_.at(currentChannel);
|
| 120 |
}
|
| 121 |
}
|
| 122 |
|
| 123 |
void
|
| 124 |
OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
|
| 125 |
{
|
| 126 |
// Retrieve necessary collections from the event.
|
| 127 |
edm::Handle<BNmuonCollection> muons;
|
| 128 |
event.getByLabel (muons_, muons);
|
| 129 |
edm::Handle<BNelectronCollection> electrons;
|
| 130 |
event.getByLabel (electrons_, electrons);
|
| 131 |
|
| 132 |
|
| 133 |
|
| 134 |
map < string, std::vector < std::vector<bool> > > electronIndividualFlags; //for each channel and each cut, vector of electrons that pass that cut
|
| 135 |
map < string, std::vector < std::vector<bool> > > electronCumulativeFlags; //for each channel and each cut, vector of electrons that pass all cuts up to and including that cut
|
| 136 |
map < string, vector <int> > passingElectronsCounter; //for each channel and each cut, number of electrons that pass all cuts up to and including that cut
|
| 137 |
|
| 138 |
map < string, std::vector < std::vector<bool> > > muonIndividualFlags; //for each channel and each cut, vector of muons that pass that cut
|
| 139 |
map < string, std::vector < std::vector<bool> > > muonCumulativeFlags; //for each channel and each cut, vector of muons that pass all cuts up to and including that cut
|
| 140 |
map < string, vector <int> > passingMuonsCounter; //for each channel and each cut, number of muons that pass all cuts up to and including that cut
|
| 141 |
|
| 142 |
|
| 143 |
//loop over all channels
|
| 144 |
for(uint currentChannel = 0; currentChannel != channels.size(); currentChannel++){
|
| 145 |
string currentChannelName = channels.at(currentChannel);
|
| 146 |
|
| 147 |
|
| 148 |
//loop over all cuts
|
| 149 |
for(uint currentCut = 0; currentCut != cutNames[currentChannelName].size(); currentCut++){
|
| 150 |
|
| 151 |
electronIndividualFlags[currentChannelName].push_back (vector<bool> ());
|
| 152 |
electronCumulativeFlags[currentChannelName].push_back (vector<bool> ());
|
| 153 |
muonIndividualFlags[currentChannelName].push_back (vector<bool> ());
|
| 154 |
muonCumulativeFlags[currentChannelName].push_back (vector<bool> ());
|
| 155 |
|
| 156 |
|
| 157 |
|
| 158 |
|
| 159 |
//loop over all electrons and set flags
|
| 160 |
for (BNelectronCollection::const_iterator electron = electrons->begin (); electron != electrons->end (); electron++){
|
| 161 |
|
| 162 |
int electronIndex = electron-electrons->begin();
|
| 163 |
bool decision = true;//electron passes if this cut doesn't cut on electrons
|
| 164 |
|
| 165 |
if(inputCollections[currentChannelName].at(currentCut) == "electrons"){
|
| 166 |
|
| 167 |
string cutVariable = cutVariables[currentChannelName].at(currentCut);
|
| 168 |
double value = valueLookup(electron, cutVariable);
|
| 169 |
|
| 170 |
decision = evaluateComparison(value,cutComparativeOperators[currentChannelName].at(currentCut),cutValues[currentChannelName].at(currentCut));
|
| 171 |
|
| 172 |
}
|
| 173 |
electronIndividualFlags[currentChannelName].at(currentCut).push_back(decision);
|
| 174 |
|
| 175 |
//set flags for electrons that pass each cut AND all the previous cuts
|
| 176 |
bool previousCumulativeFlag = true;
|
| 177 |
for(uint previousCut = 0; previousCut != currentCut; previousCut++){
|
| 178 |
if(previousCumulativeFlag && electronIndividualFlags[currentChannelName].at(previousCut).at(electronIndex)) previousCumulativeFlag = true;
|
| 179 |
else{ previousCumulativeFlag = false; break;}
|
| 180 |
}
|
| 181 |
electronCumulativeFlags[currentChannelName].at(currentCut).push_back(previousCumulativeFlag && decision);
|
| 182 |
}
|
| 183 |
|
| 184 |
|
| 185 |
//loop over all muons and set flags
|
| 186 |
for (BNmuonCollection::const_iterator muon = muons->begin (); muon != muons->end (); muon++){
|
| 187 |
|
| 188 |
int muonIndex = muon-muons->begin();
|
| 189 |
bool decision = true;//electron passes if this cut doesn't cut on electrons
|
| 190 |
|
| 191 |
//set flag to true if cut doesn't cut on muons
|
| 192 |
if(inputCollections[currentChannelName].at(currentCut) == "muons"){
|
| 193 |
|
| 194 |
string cutVariable = cutVariables[currentChannelName].at(currentCut);
|
| 195 |
double value = valueLookup(muon, cutVariable);
|
| 196 |
|
| 197 |
decision = evaluateComparison(value,cutComparativeOperators[currentChannelName].at(currentCut),cutValues[currentChannelName].at(currentCut));
|
| 198 |
|
| 199 |
}
|
| 200 |
muonIndividualFlags[currentChannelName].at(currentCut).push_back(decision);
|
| 201 |
|
| 202 |
//set flags for muons that pass each cut AND all the previous cuts
|
| 203 |
bool previousCumulativeFlag = true;
|
| 204 |
for(uint previousCut = 0; previousCut != currentCut; previousCut++){
|
| 205 |
if(previousCumulativeFlag && muonIndividualFlags[currentChannelName].at(previousCut).at(muonIndex)) previousCumulativeFlag = true;
|
| 206 |
else{ previousCumulativeFlag = false; break;}
|
| 207 |
}
|
| 208 |
muonCumulativeFlags[currentChannelName].at(currentCut).push_back(previousCumulativeFlag && decision);
|
| 209 |
}
|
| 210 |
|
| 211 |
|
| 212 |
|
| 213 |
}//end loop over all cuts
|
| 214 |
|
| 215 |
|
| 216 |
//use cumulative flags to apply cuts at event level
|
| 217 |
for(uint currentCut = 0; currentCut != cutNames[currentChannelName].size(); currentCut++){
|
| 218 |
//loop over all objects and count how many passed the cumulative selection up to this point
|
| 219 |
int numberPassing = 0;
|
| 220 |
if(inputCollections[currentChannelName].at(currentCut) == "electrons"){
|
| 221 |
for (uint electron = 0; electron != electronCumulativeFlags[currentChannelName].at(currentCut).size() ; electron++){
|
| 222 |
if(electronCumulativeFlags[currentChannelName].at(currentCut).at(electron)) numberPassing++;
|
| 223 |
}
|
| 224 |
}
|
| 225 |
else if(inputCollections[currentChannelName].at(currentCut) == "muons"){
|
| 226 |
for (uint muon = 0; muon != muonCumulativeFlags[currentChannelName].at(currentCut).size() ; muon++){
|
| 227 |
if(muonCumulativeFlags[currentChannelName].at(currentCut).at(muon)) numberPassing++;
|
| 228 |
}
|
| 229 |
}
|
| 230 |
bool cutDecision = evaluateComparison(numberPassing,eventComparativeOperators[currentChannelName].at(currentCut),numberRequiredInts[currentChannelName].at(currentCut));
|
| 231 |
cutFlows_.at(currentChannel)->at (cutNames[currentChannelName].at(currentCut)) = cutDecision;
|
| 232 |
}
|
| 233 |
|
| 234 |
cutFlows_.at(currentChannel)->fillCutFlow();
|
| 235 |
|
| 236 |
string lastCutName = cutNames[currentChannelName].at(cutNames[currentChannelName].size()-1);
|
| 237 |
if( !cutFlows_.at(currentChannel)->at (lastCutName) ) continue; //don't fill histograms if event didn't pass
|
| 238 |
|
| 239 |
vector<bool> electronFlags = electronCumulativeFlags[currentChannelName].at(cutNames[currentChannelName].size()-1);
|
| 240 |
int electronCounter = 0;
|
| 241 |
for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
|
| 242 |
if(!electronFlags.at(electronIndex)) continue;
|
| 243 |
electronCounter++;
|
| 244 |
oneDHists_.at(currentChannel)["electronPt"]->Fill (electrons->at(electronIndex).pt);
|
| 245 |
oneDHists_.at(currentChannel)["electronEta"]->Fill (electrons->at(electronIndex).eta);
|
| 246 |
oneDHists_.at(currentChannel)["electronD0"]->Fill (electrons->at(electronIndex).correctedD0Vertex);
|
| 247 |
oneDHists_.at(currentChannel)["electronAbsD0"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex));
|
| 248 |
oneDHists_.at(currentChannel)["electronD0Sig"]->Fill (electrons->at(electronIndex).correctedD0Vertex / electrons->at(electronIndex).tkD0err);
|
| 249 |
oneDHists_.at(currentChannel)["electronAbsD0Sig"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex) / electrons->at(electronIndex).tkD0err);
|
| 250 |
oneDHists_.at(currentChannel)["electronIso"]->Fill (((electrons->at(electronIndex).trackIso + electrons->at(electronIndex).caloIso) / electrons->at(electronIndex).pt));
|
| 251 |
oneDHists_.at(currentChannel)["electronFbrem"]->Fill (electrons->at(electronIndex).fbrem);
|
| 252 |
}
|
| 253 |
oneDHists_.at(currentChannel)["numElectrons"]->Fill (electronCounter);
|
| 254 |
|
| 255 |
vector<bool> muonFlags = muonCumulativeFlags[currentChannelName].at(cutNames[currentChannelName].size()-1);
|
| 256 |
int muonCounter = 0;
|
| 257 |
for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
|
| 258 |
if(!muonFlags.at(muonIndex)) continue;
|
| 259 |
muonCounter++;
|
| 260 |
oneDHists_.at(currentChannel)["muonPt"]->Fill (muons->at(muonIndex).pt);
|
| 261 |
oneDHists_.at(currentChannel)["muonEta"]->Fill (muons->at(muonIndex).eta);
|
| 262 |
oneDHists_.at(currentChannel)["muonD0"]->Fill (muons->at(muonIndex).correctedD0Vertex);
|
| 263 |
oneDHists_.at(currentChannel)["muonAbsD0"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex));
|
| 264 |
oneDHists_.at(currentChannel)["muonD0Sig"]->Fill (muons->at(muonIndex).correctedD0Vertex / muons->at(muonIndex).tkD0err);
|
| 265 |
oneDHists_.at(currentChannel)["muonAbsD0Sig"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex) / muons->at(muonIndex).tkD0err);
|
| 266 |
oneDHists_.at(currentChannel)["muonIso"]->Fill (((muons->at(muonIndex).trackIso + muons->at(muonIndex).caloIso) / muons->at(muonIndex).pt));
|
| 267 |
}
|
| 268 |
oneDHists_.at(currentChannel)["numMuons"]->Fill (muonCounter);
|
| 269 |
|
| 270 |
|
| 271 |
|
| 272 |
} //end loop over channel
|
| 273 |
|
| 274 |
masterCutFlow_->fillCutFlow();
|
| 275 |
|
| 276 |
}
|
| 277 |
|
| 278 |
|
| 279 |
bool
|
| 280 |
OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
|
| 281 |
|
| 282 |
if(comparison == ">") return testValue > cutValue;
|
| 283 |
else if(comparison == ">=") return testValue >= cutValue;
|
| 284 |
else if(comparison == "<") return testValue < cutValue;
|
| 285 |
else if(comparison == "<=") return testValue <= cutValue;
|
| 286 |
else {std::cout << "WARNING: invalid comparison operator!\n"; return false;}
|
| 287 |
|
| 288 |
}
|
| 289 |
|
| 290 |
std::vector<std::string>
|
| 291 |
OSUAnalysis::splitString (string inputString){
|
| 292 |
|
| 293 |
std::stringstream stringStream(inputString);
|
| 294 |
std::istream_iterator<std::string> begin(stringStream);
|
| 295 |
std::istream_iterator<std::string> end;
|
| 296 |
std::vector<std::string> stringVector(begin, end);
|
| 297 |
return stringVector;
|
| 298 |
|
| 299 |
}
|
| 300 |
|
| 301 |
double
|
| 302 |
OSUAnalysis::valueLookup (vector<BNelectron>::const_iterator electron, string variable){
|
| 303 |
|
| 304 |
double value = 0.0;
|
| 305 |
if(variable == "pt") value = electron->pt;
|
| 306 |
else if(variable == "eta") value = abs(electron->eta);
|
| 307 |
else if(variable == "fbrem") value = electron->fbrem;
|
| 308 |
|
| 309 |
|
| 310 |
else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
|
| 311 |
|
| 312 |
|
| 313 |
return value;
|
| 314 |
}
|
| 315 |
|
| 316 |
double
|
| 317 |
OSUAnalysis::valueLookup (vector<BNmuon>::const_iterator muon, string variable){
|
| 318 |
|
| 319 |
double value = 0.0;
|
| 320 |
if(variable == "pt") value = muon->pt;
|
| 321 |
else if(variable == "eta") value = abs(muon->eta);
|
| 322 |
|
| 323 |
else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
|
| 324 |
|
| 325 |
|
| 326 |
return value;
|
| 327 |
}
|
| 328 |
|
| 329 |
|
| 330 |
DEFINE_FWK_MODULE(OSUAnalysis);
|
