ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc
Revision: 1.111
Committed: Thu Jul 25 12:24:27 2013 UTC (11 years, 9 months ago) by wulsin
Content type: text/plain
Branch: MAIN
Changes since 1.110: +12 -0 lines
Log Message:
Create cutflow histogram in constructor, so that when running over files with 0 events a cutflow histogram with correct number of bins is still created

File Contents

# Content
1 #include "OSUT3Analysis/AnaTools/plugins/OSUAnalysis.h"
2 OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
3 // Retrieve parameters from the configuration file.
4 jets_ (cfg.getParameter<edm::InputTag> ("jets")),
5 muons_ (cfg.getParameter<edm::InputTag> ("muons")),
6 electrons_ (cfg.getParameter<edm::InputTag> ("electrons")),
7 events_ (cfg.getParameter<edm::InputTag> ("events")),
8 taus_ (cfg.getParameter<edm::InputTag> ("taus")),
9 mets_ (cfg.getParameter<edm::InputTag> ("mets")),
10 tracks_ (cfg.getParameter<edm::InputTag> ("tracks")),
11 genjets_ (cfg.getParameter<edm::InputTag> ("genjets")),
12 mcparticles_ (cfg.getParameter<edm::InputTag> ("mcparticles")),
13 stops_ (cfg.getParameter<edm::InputTag> ("stops")),
14 primaryvertexs_ (cfg.getParameter<edm::InputTag> ("primaryvertexs")),
15 bxlumis_ (cfg.getParameter<edm::InputTag> ("bxlumis")),
16 photons_ (cfg.getParameter<edm::InputTag> ("photons")),
17 superclusters_ (cfg.getParameter<edm::InputTag> ("superclusters")),
18 triggers_ (cfg.getParameter<edm::InputTag> ("triggers")),
19 trigobjs_ (cfg.getParameter<edm::InputTag> ("trigobjs")),
20 puFile_ (cfg.getParameter<string> ("puFile")),
21 deadEcalFile_ (cfg.getParameter<string> ("deadEcalFile")),
22 muonSFFile_ (cfg.getParameter<string> ("muonSFFile")),
23 dataPU_ (cfg.getParameter<string> ("dataPU")),
24 electronSFID_ (cfg.getParameter<string> ("electronSFID")),
25 muonSF_ (cfg.getParameter<string> ("muonSF")),
26 dataset_ (cfg.getParameter<string> ("dataset")),
27 datasetType_ (cfg.getParameter<string> ("datasetType")),
28 channels_ (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
29 histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
30 useEDMFormat_ (cfg.getParameter<bool>("useEDMFormat")),
31 treeBranchSets_ (cfg.getParameter<vector<edm::ParameterSet> >("treeBranchSets")),
32 plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
33 doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting")),
34 applyLeptonSF_ (cfg.getParameter<bool> ("applyLeptonSF")),
35 applyBtagSF_ (cfg.getParameter<bool> ("applyBtagSF")),
36 minBtag_ (cfg.getParameter<int> ("minBtag")),
37 printEventInfo_ (cfg.getParameter<bool> ("printEventInfo")),
38 printAllTriggers_ (cfg.getParameter<bool> ("printAllTriggers")),
39 useTrackCaloRhoCorr_ (cfg.getParameter<bool> ("useTrackCaloRhoCorr")),
40 stopCTau_ (cfg.getParameter<vector<double> > ("stopCTau")),
41 GetPlotsAfterEachCut_ (cfg.getParameter<bool> ("GetPlotsAfterEachCut")),
42 verbose_ (cfg.getParameter<int> ("verbose"))
43 {
44
45 if (verbose_) printEventInfo_ = true;
46 if (verbose_) clog << "Beginning OSUAnalysis::OSUAnalysis constructor." << endl;
47
48 TH1::SetDefaultSumw2();
49
50 //create pile-up reweighting object, if necessary
51 if(datasetType_ != "data") {
52 if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
53 if (applyLeptonSF_){
54 muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
55 electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
56 }
57 if (applyBtagSF_){
58 bTagSFWeight_ = new BtagSFWeight;
59 }
60 }
61 if (datasetType_ == "signalMC" && regex_match (dataset_, regex ("stop.*to.*_.*mm.*")))
62 stopCTauWeight_ = new StopCTauWeight (stopCTau_.at(0), stopCTau_.at(1), stops_);
63
64
65 // Construct Cutflow Objects. These store the results of cut decisions and
66 // handle filling cut flow histograms.
67 masterCutFlow_ = new CutFlow (fs_);
68
69 //always get vertex collection so we can assign the primary vertex in the event
70 objectsToGet.push_back("primaryvertexs");
71 objectsToPlot.push_back("primaryvertexs");
72 objectsToFlag.push_back("primaryvertexs");
73
74
75 //always get the MC particles to do GEN-matching
76 objectsToGet.push_back("mcparticles");
77
78 //always get the event collection to do pile-up reweighting
79 objectsToGet.push_back("events");
80
81
82 // Parse the tree variable definitions.
83 for (uint iBranchSet = 0; !useEDMFormat_ && iBranchSet<treeBranchSets_.size(); iBranchSet++) {
84 string tempInputCollection = treeBranchSets_.at(iBranchSet).getParameter<string> ("inputCollection");
85 if(tempInputCollection.find("pairs")!=string::npos) { clog << "Warning: tree filling is not configured for pairs of objects, so will not work for collection: " << tempInputCollection << endl; }
86 objectsToGet.push_back(tempInputCollection);
87 objectsToFlag.push_back(tempInputCollection);
88
89 vector<string> branchList(treeBranchSets_.at(iBranchSet).getParameter<vector<string> >("branches"));
90
91 for (uint iBranch = 0; iBranch<branchList.size(); iBranch++) {
92 BranchSpecs br;
93 br.inputCollection = tempInputCollection;
94 br.inputVariable = branchList.at(iBranch);
95 TString newName = TString(br.inputCollection) + "_" + TString(br.inputVariable);
96 br.name = string(newName.Data());
97 treeBranches_.push_back(br);
98 if (verbose_>3) clog << " Adding branch to BNTree: " << br.name << endl;
99 }
100
101 } // end for (uint iBranchSet = 0; iBranchSet<treeBranchSets_.size(); iBranchSet++)
102
103
104 //parse the histogram definitions
105 for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
106
107 string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
108 if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
109 if(tempInputCollection == "photon-muon pairs") tempInputCollection = "muon-photon pairs";
110 if(tempInputCollection == "photon-electron pairs") tempInputCollection = "electron-photon pairs";
111 if(tempInputCollection == "jet-electron pairs") tempInputCollection = "electron-jet pairs";
112 if(tempInputCollection == "jet-photon pairs") tempInputCollection = "photon-jet pairs";
113 if(tempInputCollection == "jet-muon pairs") tempInputCollection = "muon-jet pairs";
114 if(tempInputCollection == "event-muon pairs") tempInputCollection = "muon-event pairs";
115 if(tempInputCollection == "jet-met pairs") tempInputCollection = "met-jet pairs";
116 if(tempInputCollection == "track-jet pairs") tempInputCollection = "track-jet pairs";
117 if(tempInputCollection == "event-track pairs") tempInputCollection = "track-event pairs";
118 if(tempInputCollection == "secondary muon-muon pairs") tempInputCollection = "muon-secondary muon pairs";
119 if(tempInputCollection == "secondary jet-muon pairs") tempInputCollection = "muon-secondary jet pairs";
120 if(tempInputCollection == "secondary photon-muon pairs") tempInputCollection = "muon-secondary photon pairs";
121 if(tempInputCollection == "secondary jet-electron pairs") tempInputCollection = "electron-secondary jet pairs";
122 if(tempInputCollection == "secondary jet-photon pairs") tempInputCollection = "photon-secondary jet pairs";
123 if(tempInputCollection == "secondary jet-jet pairs") tempInputCollection = "jet-secondary jet pairs";
124 if(tempInputCollection == "secondary electron-electron pairs") tempInputCollection = "electron-secondary electron pairs";
125 if(tempInputCollection == "trigobj-electron pairs") tempInputCollection = "electron-trigobj pairs";
126 if(tempInputCollection == "trigobj-muon pairs") tempInputCollection = "muon-trigobj pairs";
127 if(tempInputCollection.find("pairs")==string::npos){ //just a single object
128 if(tempInputCollection.find("secondary")!=string::npos){//secondary object
129 int spaceIndex = tempInputCollection.find(" ");
130 int secondWordLength = tempInputCollection.size() - spaceIndex;
131 objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
132 }
133 else{
134 objectsToGet.push_back(tempInputCollection);
135 }
136 objectsToPlot.push_back(tempInputCollection);
137 objectsToFlag.push_back(tempInputCollection);
138 } else { //pair of objects, need to add the pair and the individual objects to the lists of things to Get/Plot/Cut
139 string obj1;
140 string obj2;
141 getTwoObjs(tempInputCollection, obj1, obj2);
142 string obj2ToGet = getObjToGet(obj2);
143 objectsToFlag.push_back(tempInputCollection);
144 objectsToFlag.push_back(obj1);
145 objectsToFlag.push_back(obj2);
146 objectsToPlot.push_back(tempInputCollection);
147 objectsToPlot.push_back(obj1);
148 objectsToPlot.push_back(obj2);
149 objectsToGet.push_back(tempInputCollection);
150 objectsToGet.push_back(obj1);
151 objectsToGet.push_back(obj2ToGet);
152 } // end else
153 if (find(objectsToPlot.begin(), objectsToPlot.end(), "events") != objectsToPlot.end())
154 objectsToGet.push_back("jets");
155
156
157 vector<edm::ParameterSet> histogramList_ (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
158
159 for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
160
161 vector<double> defaultValue;
162 defaultValue.push_back (-1.0);
163
164 histogram tempHistogram;
165 tempHistogram.inputCollection = tempInputCollection;
166 tempHistogram.name = histogramList_.at(currentHistogram).getParameter<string>("name");
167 tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
168 tempHistogram.bins = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("bins", defaultValue);
169 tempHistogram.variableBinsX = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("variableBinsX", defaultValue);
170 tempHistogram.variableBinsY = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("variableBinsY", defaultValue);
171 tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
172
173 histograms.push_back(tempHistogram);
174
175 }
176 } // for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++)
177
178 //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
179 sort( objectsToPlot.begin(), objectsToPlot.end() );
180 objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
181
182
183
184 //add histograms with the gen-matched id, mother id, and grandmother id
185 for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
186
187 string currentObject = objectsToPlot.at(currentObjectIndex);
188 if(currentObject != "muons" && currentObject != "secondary muons" && currentObject != "secondary electrons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "secondary photons"&& currentObject != "superclusters") continue;
189
190 histogram tempIdHisto;
191 histogram tempMomIdHisto;
192 histogram tempGmaIdHisto;
193 histogram tempIdVsMomIdHisto;
194 histogram tempIdVsGmaIdHisto;
195
196 tempIdHisto.inputCollection = currentObject;
197 tempMomIdHisto.inputCollection = currentObject;
198 tempGmaIdHisto.inputCollection = currentObject;
199 tempIdVsMomIdHisto.inputCollection = currentObject;
200 tempIdVsGmaIdHisto.inputCollection = currentObject;
201
202 if(currentObject == "secondary muons") currentObject = "secondaryMuons";
203 if(currentObject == "secondary photons") currentObject = "secondaryPhotons";
204 if(currentObject == "secondary electrons") currentObject = "secondaryElectrons";
205
206 currentObject = currentObject.substr(0, currentObject.size()-1);
207 tempIdHisto.name = currentObject+"GenMatchId";
208 tempMomIdHisto.name = currentObject+"GenMatchMotherId";
209 tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
210 tempIdVsMomIdHisto.name = currentObject+"GenMatchIdVsMotherId";
211 tempIdVsGmaIdHisto.name = currentObject+"GenMatchIdVsGrandmotherId";
212
213 currentObject.at(0) = toupper(currentObject.at(0));
214 tempIdHisto.title = currentObject+" Gen-matched Particle";
215 tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
216 tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
217 tempIdVsMomIdHisto.title = currentObject+" Gen-matched Particle's Mother vs. Particle;Particle;Mother";
218 tempIdVsGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother vs. Particle;Particle;Grandmother";
219
220
221 int maxNum = 25;
222 vector<double> binVector;
223 binVector.push_back(maxNum);
224 binVector.push_back(0);
225 binVector.push_back(maxNum);
226
227 tempIdHisto.bins = binVector;
228 tempIdHisto.inputVariables.push_back("genMatchedId");
229 tempMomIdHisto.bins = binVector;
230 tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
231 tempGmaIdHisto.bins = binVector;
232 tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
233 binVector.push_back(maxNum);
234 binVector.push_back(0);
235 binVector.push_back(maxNum);
236 tempIdVsMomIdHisto.bins = binVector;
237 tempIdVsMomIdHisto.inputVariables.push_back("genMatchedId");
238 tempIdVsMomIdHisto.inputVariables.push_back("genMatchedMotherIdReverse");
239 tempIdVsGmaIdHisto.bins = binVector;
240 tempIdVsGmaIdHisto.inputVariables.push_back("genMatchedId");
241 tempIdVsGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherIdReverse");
242
243 histograms.push_back(tempIdHisto);
244 histograms.push_back(tempMomIdHisto);
245 histograms.push_back(tempGmaIdHisto);
246 histograms.push_back(tempIdVsMomIdHisto);
247 histograms.push_back(tempIdVsGmaIdHisto);
248 }
249
250
251 channel tempChannel;
252 //loop over all channels (event selections)
253 for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
254
255 //get name of channel
256 string channelName (channels_.at(currentChannel).getParameter<string>("name"));
257 tempChannel.name = channelName;
258 TString channelLabel = channelName;
259 if (verbose_) clog << "Processing channel: " << channelName << endl;
260
261 //set triggers for this channel
262 vector<string> triggerNames;
263 triggerNames.clear();
264 vector<string> triggerToVetoNames;
265 triggerToVetoNames.clear();
266
267 tempChannel.triggers.clear();
268 tempChannel.triggersToVeto.clear();
269 if(channels_.at(currentChannel).exists("triggers")){
270 triggerNames = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
271 for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
272 tempChannel.triggers.push_back(triggerNames.at(trigger));
273 objectsToGet.push_back("triggers");
274 }
275 if(channels_.at(currentChannel).exists("triggersToVeto")){
276 triggerToVetoNames = channels_.at(currentChannel).getParameter<vector<string> >("triggersToVeto");
277 for(uint trigger = 0; trigger!= triggerToVetoNames.size(); trigger++)
278 tempChannel.triggersToVeto.push_back(triggerToVetoNames.at(trigger));
279 objectsToGet.push_back("triggers");
280 }
281
282
283
284 //create cutFlow for this channel
285 cutFlows_.push_back (new CutFlow (fs_, channelName));
286 vector<TFileDirectory> directories; //vector of directories in the output file.
287 vector<TFileDirectory> treeDirectories; //vector of directories for trees in the output file.
288 vector<string> subSubDirNames;//subdirectories in each channel.
289 //get list of cuts for this channel
290 vector<edm::ParameterSet> cuts_ (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
291
292
293 //loop over and parse all cuts
294 for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
295 cut tempCut;
296 //store input collection for cut
297 string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
298 if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
299 if(tempInputCollection == "photon-electron pairs") tempInputCollection = "electron-photon pairs";
300 if(tempInputCollection == "photon-muon pairs") tempInputCollection = "muon-photon pairs";
301 if(tempInputCollection == "jet-electron pairs") tempInputCollection = "electron-jet pairs";
302 if(tempInputCollection == "jet-muon pairs") tempInputCollection = "muon-jet pairs";
303 if(tempInputCollection == "event-muon pairs") tempInputCollection = "muon-event pairs";
304 if(tempInputCollection == "jet-met pairs") tempInputCollection = "met-jet pairs";
305 if(tempInputCollection == "track-jet pairs") tempInputCollection = "track-jet pairs";
306 if(tempInputCollection == "jet-photon pairs") tempInputCollection = "photon-jet pairs";
307 if(tempInputCollection == "event-track pairs") tempInputCollection = "track-event pairs";
308 if(tempInputCollection == "secondary muon-muon pairs") tempInputCollection = "muon-secondary muon pairs";
309 if(tempInputCollection == "secondary jet-jet pairs") tempInputCollection = "jet-secondary jet pairs";
310 if(tempInputCollection == "secondary jet-muon pairs") tempInputCollection = "muon-secondary jet pairs";
311 if(tempInputCollection == "secondary photon-muon pairs") tempInputCollection = "muon-secondary photon pairs";
312 if(tempInputCollection == "secondary jet-photon pairs") tempInputCollection = "photon-secondary jet pairs";
313 if(tempInputCollection == "secondary jet-electron pairs") tempInputCollection = "electron-secondary jet pairs";
314 if(tempInputCollection == "secondary electron-electron pairs") tempInputCollection = "electron-secondary electron pairs";
315 if(tempInputCollection == "trigobj-electron pairs") tempInputCollection = "electron-trigobj pairs";
316 if(tempInputCollection == "trigobj-muon pairs") tempInputCollection = "muon-trigobj pairs";
317 tempCut.inputCollection = tempInputCollection;
318 if(tempInputCollection.find("pairs")==string::npos){ //just a single object
319 if(tempInputCollection.find("secondary")!=string::npos){//secondary object
320 int spaceIndex = tempInputCollection.find(" ");
321 int secondWordLength = tempInputCollection.size() - spaceIndex;
322 objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
323 }
324 else{
325 objectsToGet.push_back(tempInputCollection);
326 }
327 objectsToCut.push_back(tempInputCollection);
328 objectsToFlag.push_back(tempInputCollection);
329 }
330 else{//pair of objects, need to add them both to objectsToGet
331 string obj1;
332 string obj2;
333 getTwoObjs(tempInputCollection, obj1, obj2);
334 string obj2ToGet = getObjToGet(obj2);
335 objectsToCut.push_back(tempInputCollection);
336 objectsToCut.push_back(obj1);
337 objectsToCut.push_back(obj2);
338 objectsToFlag.push_back(tempInputCollection);
339 objectsToFlag.push_back(obj1);
340 objectsToFlag.push_back(obj2);
341 objectsToGet.push_back(tempInputCollection);
342 objectsToGet.push_back(obj1);
343 objectsToGet.push_back(obj2ToGet);
344
345 }
346
347
348
349 //split cut string into parts and store them
350 string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
351 vector<string> cutStringVector = splitString(cutString);
352 if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
353 clog << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
354 exit(0);
355 }
356 tempCut.numSubcuts = (cutStringVector.size()+1)/4;
357 for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
358 int indexOffset = 4 * subcutIndex;
359 string currentVariableString = cutStringVector.at(indexOffset);
360 if(currentVariableString.find("(")==string::npos){
361 tempCut.functions.push_back("");//no function was specified
362 tempCut.variables.push_back(currentVariableString);// variable to cut on
363 }
364 else{
365 tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
366 string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
367 tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
368 }
369 tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
370 tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
371 tempCut.cutStringValues.push_back(cutStringVector.at(indexOffset+2));// string value to pass cut
372 if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
373 }
374
375 //get number of objects required to pass cut for event to pass
376 string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
377 vector<string> numberRequiredVector = splitString(numberRequiredString);
378 if(numberRequiredVector.size()!=2){
379 clog << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
380 exit(0);
381 }
382
383 int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
384 tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
385 tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
386
387 //Set up vectors to store the directories and subDirectories for each channel.
388 string subSubDirName;
389 string tempCutName;
390 if(cuts_.at(currentCut).exists("alias")){
391 tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
392 subSubDirName = "After " + tempCutName + " Cut Applied";
393 }
394 else{
395 //construct string for cutflow table
396 bool plural = numberRequiredInt != 1;
397 string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
398 string cutName = numberRequiredString + " " + collectionString + " with " + cutString;
399 tempCutName = cutName;
400 subSubDirName = "After " + numberRequiredString + " " + collectionString + " with " + cutString + " Cut Applied";
401 }
402
403 tempCut.isVeto = false;
404 if(cuts_.at(currentCut).exists("isVeto")){
405 bool isVeto = cuts_.at(currentCut).getParameter<bool> ("isVeto");
406 if (isVeto == true) tempCut.isVeto = true;
407 }
408 subSubDirNames.push_back(subSubDirName);
409 tempCut.name = tempCutName;
410 if (verbose_) clog << "Setting up cut, index: " << currentCut << ", input collection: " << tempInputCollection<< ", name: " << tempCut.name << endl;
411
412 tempChannel.cuts.push_back(tempCut);
413
414
415 }//end loop over cuts
416
417 vector<map<string, TH1D*>> oneDHistsTmp;
418 vector<map<string, TH2D*>> twoDHistsTmp;
419 //book a directory in the output file with the name of the channel
420 TFileDirectory subDir = fs_->mkdir( channelName );
421 //loop over the cuts to set up subdirectory for each cut if GetPlotsAfterEachCut_ is true.
422 if(GetPlotsAfterEachCut_){
423 for( uint currentDir=0; currentDir != subSubDirNames.size(); currentDir++){
424 TFileDirectory subSubDir = subDir.mkdir( subSubDirNames[currentDir] );
425 directories.push_back(subSubDir);
426 map<string, TH1D*> oneDhistoMap;
427 oneDHistsTmp.push_back(oneDhistoMap);
428 map<string, TH2D*> twoDhistoMap;
429 twoDHistsTmp.push_back(twoDhistoMap);
430 }
431 treeDirectories.push_back(subDir);
432 oneDHists_.push_back(oneDHistsTmp);
433 twoDHists_.push_back(twoDHistsTmp);
434 }
435 //only set up directories with names of the channels if GetPlotsAfterEachCut_ is false.
436 else{
437 map<string, TH1D*> oneDhistoMap;
438 oneDHistsTmp.push_back(oneDhistoMap);
439 map<string, TH2D*> twoDhistoMap;
440 twoDHistsTmp.push_back(twoDhistoMap);
441 oneDHists_.push_back(oneDHistsTmp);
442 twoDHists_.push_back(twoDHistsTmp);
443 directories.push_back(subDir);
444 treeDirectories.push_back(subDir);
445
446 }
447
448 for(uint currentDir = 0; !useEDMFormat_ && currentDir != treeDirectories.size(); currentDir++){
449 TTree* newTree = treeDirectories.at(currentDir).make<TTree> (TString("BNTree_"+channelLabel), TString("BNTree_"+channelLabel));
450 BNTrees_.push_back(newTree);
451 BNTrees_.back()->Branch("event_runInt", &BNTreeBranchVals_runInt_, "event_runInt/I");
452 BNTrees_.back()->Branch("event_lumiInt", &BNTreeBranchVals_lumiInt_, "event_lumiInt/I");
453 BNTrees_.back()->Branch("event_evtLong", &BNTreeBranchVals_evtLong_, "event_evtLong/L");
454 for (uint iBranch = 0; iBranch < treeBranches_.size(); iBranch++){
455 BranchSpecs currentVar = treeBranches_.at(iBranch);
456 vector<float> newVec;
457 BNTreeBranchVals_[currentVar.name] = newVec;
458 BNTrees_.back()->Branch(TString(currentVar.name), &BNTreeBranchVals_.at(currentVar.name));
459 } // end for (uint iBranch = 0; iBranch < treeBranches_.size(); iBranch++)
460 }
461
462 //book all histograms included in the configuration
463 for(uint currentDir = 0; currentDir != directories.size(); currentDir++){//loop over all the directories.
464
465 for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
466
467 histogram currentHistogram = histograms.at(currentHistogramIndex);
468 int numBinsElements = currentHistogram.bins.size();
469 int numInputVariables = currentHistogram.inputVariables.size();
470 int numBinEdgesX = currentHistogram.variableBinsX.size();
471 int numBinEdgesY = currentHistogram.variableBinsY.size();
472
473 if(numBinsElements == 1){
474 if(numBinEdgesX > 1){
475 if(numBinEdgesY > 1)
476 numBinsElements = 6;
477 else
478 numBinsElements = 3;
479 }
480 }
481 if(numBinsElements != 3 && numBinsElements !=6){
482 clog << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
483 exit(0);
484 }
485 else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
486 clog << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
487 exit(0);
488 }
489 else if(numBinsElements == 3){
490 if (currentHistogram.bins.size () == 3)
491 oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH1D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, currentHistogram.bins.at(0), currentHistogram.bins.at(1), currentHistogram.bins.at(2));
492 else
493 {
494 oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH1D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, numBinEdgesX - 1, currentHistogram.variableBinsX.data ());
495 }
496 }
497 else if(numBinsElements == 6){
498 if (currentHistogram.bins.size () == 6)
499 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH2D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, currentHistogram.bins.at(0), currentHistogram.bins.at(1), currentHistogram.bins.at(2),currentHistogram.bins.at(3),currentHistogram.bins.at(4),currentHistogram.bins.at(5));
500 else
501 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH2D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, numBinEdgesX - 1, currentHistogram.variableBinsX.data (), numBinEdgesY - 1, currentHistogram.variableBinsY.data ());
502 }
503
504
505 if(currentHistogram.name.find("GenMatch")==string::npos) continue;
506
507 // bin particle type
508 // --- -------------
509 // 0 unmatched
510 // 1 u
511 // 2 d
512 // 3 s
513 // 4 c
514 // 5 b
515 // 6 t
516 // 7 e
517 // 8 mu
518 // 9 tau
519 // 10 nu
520 // 11 g
521 // 12 gamma
522 // 13 Z
523 // 14 W
524 // 15 light meson
525 // 16 K meson
526 // 17 D meson
527 // 18 B meson
528 // 19 light baryon
529 // 20 strange baryon
530 // 21 charm baryon
531 // 22 bottom baryon
532 // 23 QCD string
533 // 24 other
534
535 vector<TString> labelArray;
536 labelArray.push_back("unmatched");
537 labelArray.push_back("u");
538 labelArray.push_back("d");
539 labelArray.push_back("s");
540 labelArray.push_back("c");
541 labelArray.push_back("b");
542 labelArray.push_back("t");
543 labelArray.push_back("e");
544 labelArray.push_back("#mu");
545 labelArray.push_back("#tau");
546 labelArray.push_back("#nu");
547 labelArray.push_back("g");
548 labelArray.push_back("#gamma");
549 labelArray.push_back("Z");
550 labelArray.push_back("W");
551 labelArray.push_back("light meson");
552 labelArray.push_back("K meson");
553 labelArray.push_back("D meson");
554 labelArray.push_back("B meson");
555 labelArray.push_back("light baryon");
556 labelArray.push_back("strange baryon");
557 labelArray.push_back("charm baryon");
558 labelArray.push_back("bottom baryon");
559 labelArray.push_back("QCD string");
560 labelArray.push_back("other");
561
562 for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
563 if(currentHistogram.name.find("GenMatchIdVsMotherId")==string::npos && currentHistogram.name.find("GenMatchIdVsGrandmotherId")==string::npos) {
564 oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
565 }
566 else {
567 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetYaxis()->SetBinLabel(bin+1,labelArray.at(currentHistogram.bins.at(0)-bin-1));
568 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
569 }
570 }
571 if(currentHistogram.name.find("GenMatchIdVsMotherId")!=string::npos || currentHistogram.name.find("GenMatchIdVsGrandmotherId")!=string::npos) {
572 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->CenterTitle();
573 twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetYaxis()->CenterTitle();
574 }
575
576 } // end for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++)
577
578
579 // Book a histogram for the number of each object type to be plotted.
580 // Name of objectToPlot here must match the name specified in OSUAnalysis::analyze().
581 for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
582 string currentObject = objectsToPlot.at(currentObjectIndex);
583 int maxNum = 10;
584 if(currentObject == "mcparticles") maxNum = 50;
585 else if(currentObject == "primaryvertexs") maxNum = 50;
586
587 if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
588 else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
589 else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
590 else if(currentObject == "electron-photon pairs") currentObject = "electronPhotonPairs";
591 else if(currentObject == "muon-photon pairs") currentObject = "muonPhotonPairs";
592 else if(currentObject == "secondary jets") currentObject = "secondaryJets";
593 else if(currentObject == "secondary photons") currentObject = "secondaryPhotons";
594 else if(currentObject == "jet-jet pairs") currentObject = "dijetPairs";
595 else if(currentObject == "jet-secondary jet pairs") currentObject = "jetSecondaryJetPairs";
596 else if(currentObject == "electron-jet pairs") currentObject = "electronJetPairs";
597 else if(currentObject == "muon-jet pairs") currentObject = "muonJetPairs";
598 else if(currentObject == "muon-event pairs") currentObject = "muonEventPairs";
599 else if(currentObject == "photon-jet pairs") currentObject = "photonJetPairs";
600 else if(currentObject == "met-jet pairs") currentObject = "metJetPairs";
601 else if(currentObject == "track-jet pairs") currentObject = "trackJetPairs";
602 else if(currentObject == "muon-secondary jet pairs") currentObject = "muonSecondaryJetPairs";
603 else if(currentObject == "muon-secondary photon pairs") currentObject = "muonSecondaryPhotonPairs";
604 else if(currentObject == "photon-secondary jet pairs") currentObject = "photonSecondaryJetPairs";
605 else if(currentObject == "electron-secondary jet pairs") currentObject = "electronSecondaryJetPairs";
606 else if(currentObject == "track-event pairs") currentObject = "trackEventPairs";
607 else if(currentObject == "electron-track pairs") currentObject = "electronTrackPairs";
608 else if(currentObject == "muon-track pairs") currentObject = "muonTrackPairs";
609 else if(currentObject == "muon-tau pairs") currentObject = "muonTauPairs";
610 else if(currentObject == "tau-tau pairs") currentObject = "ditauPairs";
611 else if(currentObject == "tau-track pairs") currentObject = "tauTrackPairs";
612 else if(currentObject == "muon-secondary muon pairs") currentObject = "muonSecondaryMuonPairs";
613 else if(currentObject == "secondary muons") currentObject = "secondaryMuons";
614 else if(currentObject == "electron-secondary electron pairs") currentObject = "electronSecondaryElectronPairs";
615 else if(currentObject == "secondary electrons") currentObject = "secondaryElectrons";
616 else if(currentObject == "electron-trigobj pairs") currentObject = "electronTrigobjPairs";
617 else if(currentObject == "muon-trigobj pairs") currentObject = "muonTrigobjPairs";
618
619 currentObject.at(0) = toupper(currentObject.at(0));
620 string histoName = "num" + currentObject;
621
622 if(histoName == "numPrimaryvertexs"){
623 string newHistoName = histoName + "BeforePileupCorrection";
624 oneDHists_.at(currentChannel).at(currentDir)[newHistoName] = directories.at(currentDir).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
625 newHistoName = histoName + "AfterPileupCorrection";
626 oneDHists_.at(currentChannel).at(currentDir)[newHistoName] = directories.at(currentDir).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
627 }
628 else
629 oneDHists_.at(currentChannel).at(currentDir)[histoName] = directories.at(currentDir).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
630 }
631
632 if (verbose_) {
633 clog << "List of 1D histograms booked for channel " << channelName << " and directory " << currentDir << ":" << endl;
634 for(map<string, TH1D*>::iterator
635 iter = oneDHists_.at(currentChannel).at(currentDir).begin();
636 iter != oneDHists_.at(currentChannel).at(currentDir).end();
637 iter++) {
638 clog << " " << iter->first
639 << ": name=" << iter->second->GetName()
640 << ", title=" << iter->second->GetTitle()
641 << ", bins=(" << iter->second->GetNbinsX()
642 << "," << iter->second->GetXaxis()->GetXmin()
643 << "," << iter->second->GetXaxis()->GetXmax()
644 << ")" << endl;
645 }
646 clog << "List of 2D histograms booked for channel " << channelName << " and directory " << currentDir << ":" << endl;
647 for(map<string, TH2D*>::iterator
648 iter = twoDHists_.at(currentChannel).at(currentDir).begin();
649 iter != twoDHists_.at(currentChannel).at(currentDir).end();
650 iter++) {
651 clog << " " << iter->first
652 << ": name=" << iter->second->GetName()
653 << ", title=" << iter->second->GetTitle()
654 << ", binsX=(" << iter->second->GetNbinsX()
655 << "," << iter->second->GetXaxis()->GetXmin()
656 << "," << iter->second->GetXaxis()->GetXmax()
657 << ")"
658 << ", binsY=(" << iter->second->GetNbinsY()
659 << "," << iter->second->GetYaxis()->GetXmin()
660 << "," << iter->second->GetYaxis()->GetXmax()
661 << ")"
662 << endl;
663 }
664 }
665
666 }//end of loop over directories
667 channels.push_back(tempChannel);
668 tempChannel.cuts.clear();
669 }//end loop over channels
670
671
672 // Create the cutflow histogram and fill with 0 weight, in case no events are found in the input file.
673 for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
674 channel currentChannel = channels.at(currentChannelIndex);
675 if(currentChannel.triggers.size() != 0 || currentChannel.triggersToVeto.size() != 0){ //triggers specified
676 cutFlows_.at(currentChannelIndex)->at("trigger") = true;
677 }
678 for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
679 cut currentCut = currentChannel.cuts.at(currentCutIndex);
680 cutFlows_.at(currentChannelIndex)->at (currentCut.name) = true;
681 }
682 cutFlows_.at(currentChannelIndex)->fillCutFlow(0); // fill cutflow with 0 weight, just to create the cutflow histograms
683 }
684
685 //make unique vector of objects we need to get from the event
686 sort( objectsToGet.begin(), objectsToGet.end() );
687 objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
688 //make unique vector of objects we need to cut on
689 sort( objectsToCut.begin(), objectsToCut.end() );
690 objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
691 //make unique vector of objects we need to set flags for
692 sort( objectsToFlag.begin(), objectsToFlag.end() );
693 objectsToFlag.erase( unique( objectsToFlag.begin(), objectsToFlag.end() ), objectsToFlag.end() );
694
695 // Move all paired objects to the end of the list, so that the flags for single objects are always set before those of paired objects.
696 for (uint i=0; i<objectsToFlag.size(); i++) {
697 if (objectsToFlag.at(i).find("pairs")!=string::npos) { // if it contains "pairs"
698 objectsToFlag.push_back(objectsToFlag.at(i));
699 objectsToFlag.erase(objectsToFlag.begin()+i);
700 }
701 }
702
703 if (verbose_) {
704 clog << "List of channels:" << endl;
705 for (uint i=0; i<channels.size(); i++) {
706 clog << " " << channels.at(i).name << endl;
707 }
708 clog << "List of objects to get:" << endl;
709 for (uint i=0; i<objectsToGet.size(); i++) {
710 clog << " " << objectsToGet.at(i) << endl;
711 }
712 clog << "List of objects to plot:" << endl;
713 for (uint i=0; i<objectsToPlot.size(); i++) {
714 clog << " " << objectsToPlot.at(i) << endl;
715 }
716 clog << "List of objects to cut:" << endl;
717 for (uint i=0; i<objectsToCut.size(); i++) {
718 clog << " " << objectsToCut.at(i) << endl;
719 }
720 clog << "List of objects to flag:" << endl;
721 for (uint i=0; i<objectsToFlag.size(); i++) {
722 clog << " " << objectsToFlag.at(i) << endl;
723 }
724 }
725
726 produces<map<string, bool> > ("channelDecisions");
727
728 if (printEventInfo_) {
729 findEventsLog = new ofstream();
730 findEventsLog->open("findEvents.txt");
731 clog << "Listing run:lumi:event in file findEvents.txt for events that pass full selection (of any channel)." << endl;
732 } else {
733 findEventsLog = 0;
734 }
735
736 isFirstEvent_ = true;
737
738 if (verbose_) clog << "Finished OSUAnalysis::OSUAnalysis constructor." << endl;
739
740
741 } // end constructor OSUAnalysis::OSUAnalysis()
742
743
744 OSUAnalysis::~OSUAnalysis ()
745 {
746
747 if (verbose_) clog << "Beginning OSUAnalysis::OSUAnalysis destructor." << endl;
748
749 // Destroying the CutFlow objects causes the cut flow numbers and time
750 // information to be printed to standard output.
751 for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
752 delete cutFlows_.at(currentChannel);
753 }
754
755 if (printEventInfo_ && findEventsLog) findEventsLog->close();
756
757 clog << "=============================================" << endl;
758 clog << "Successfully completed OSUAnalysis." << endl;
759 clog << "=============================================" << endl;
760
761 if (verbose_) clog << "Finished OSUAnalysis::OSUAnalysis destructor." << endl;
762
763 }
764
765 void
766 OSUAnalysis::produce (edm::Event &event, const edm::EventSetup &setup)
767 {
768 // Retrieve necessary collections from the event.
769
770 if (verbose_) clog << "Beginning OSUAnalysis::produce." << endl;
771
772 if (find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
773 event.getByLabel (triggers_, triggers);
774
775 if (find(objectsToGet.begin(), objectsToGet.end(), "trigobjs") != objectsToGet.end())
776 event.getByLabel (trigobjs_, trigobjs);
777
778 if (find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
779 event.getByLabel (jets_, jets);
780
781 if (find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
782 event.getByLabel (muons_, muons);
783
784 if (find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
785 event.getByLabel (electrons_, electrons);
786
787 if (find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
788 event.getByLabel (events_, events);
789
790 if (find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
791 event.getByLabel (taus_, taus);
792
793 if (find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
794 event.getByLabel (mets_, mets);
795
796 if (find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
797 event.getByLabel (tracks_, tracks);
798
799 if (find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
800 event.getByLabel (genjets_, genjets);
801
802 if (find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
803 event.getByLabel (mcparticles_, mcparticles);
804
805 if (find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
806 event.getByLabel (primaryvertexs_, primaryvertexs);
807
808 if (find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
809 event.getByLabel (bxlumis_, bxlumis);
810
811 if (find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
812 event.getByLabel (photons_, photons);
813
814 if (find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
815 event.getByLabel (superclusters_, superclusters);
816
817 if (datasetType_ == "signalMC"){
818 if (find(objectsToGet.begin(), objectsToGet.end(), "stops") != objectsToGet.end())
819 event.getByLabel (stops_, stops);
820 }
821
822 if (useTrackCaloRhoCorr_) {
823 // Used only for pile-up correction of by-hand calculation of isolation energy.
824 // This rho collection is not available in all BEANs.
825 // For description of rho values for different jet reconstruction algorithms, see
826 // https://twiki.cern.ch/twiki/bin/view/CMS/JetAlgorithms#Algorithms
827 event.getByLabel ("kt6CaloJets","rho", rhokt6CaloJetsHandle_);
828 }
829
830 double masterScaleFactor = 1.0;
831
832 //get pile-up event weight
833 if (doPileupReweighting_ && datasetType_ != "data") {
834 //for "data" datasets, the numTruePV is always set to -1
835 if (events->at(0).numTruePV < 0 && isFirstEvent_) clog << "WARNING[OSUAnalysis::analyze]: Event has numTruePV<0. Turning off pile-up reweighting." << endl;
836 else masterScaleFactor *= puWeight_->at (events->at (0).numTruePV);
837 }
838
839 stopCTauScaleFactor_ = 1.0;
840 if (datasetType_ == "signalMC" && regex_match (dataset_, regex ("stop.*to.*_.*mm.*")))
841 stopCTauScaleFactor_ = stopCTauWeight_->at (event);
842 masterScaleFactor *= stopCTauScaleFactor_;
843
844 //loop over all channels
845
846 auto_ptr<map<string, bool> > channelDecisions (new map<string, bool>);
847 for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
848 channel currentChannel = channels.at(currentChannelIndex);
849 if (verbose_>1) clog << " Processing channel " << currentChannel.name << endl;
850
851 flagMap individualFlags;
852 counterMap passingCounter;
853 cumulativeFlags.clear ();
854
855 for (map<string, vector<float>>::iterator iter = BNTreeBranchVals_.begin();
856 iter != BNTreeBranchVals_.end(); iter++) {
857 iter->second.clear(); // clear array
858 }
859
860 bool triggerDecision = true;
861 if(currentChannel.triggers.size() != 0 || currentChannel.triggersToVeto.size() != 0){ //triggers specified
862 triggerDecision = evaluateTriggers(currentChannel.triggers, currentChannel.triggersToVeto, triggers.product());
863 cutFlows_.at(currentChannelIndex)->at ("trigger") = triggerDecision;
864 }
865
866 //loop over all cuts
867 for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
868 cut currentCut = currentChannel.cuts.at(currentCutIndex);
869 if (verbose_>2) clog << " Processing cut, index: " << currentCutIndex << ", input collection: " << currentCut.inputCollection << ", name: " << currentCut.name << endl;
870
871 for(uint currentObjectIndex = 0; currentObjectIndex != objectsToFlag.size(); currentObjectIndex++){
872 string currentObject = objectsToFlag.at(currentObjectIndex);
873
874 //initialize maps to get ready to set cuts
875 individualFlags[currentObject].push_back (vector<pair<bool,bool> > ());
876 cumulativeFlags[currentObject].push_back (vector<pair<bool,bool> > ());
877
878 }
879
880 //set flags for all relevant objects
881 for(int currentObjectIndex = -1; currentObjectIndex != int(objectsToFlag.size()); currentObjectIndex++){
882 string currentObject;
883 if (currentObjectIndex < 0) currentObject = currentCut.inputCollection; // In the first loop, set the flags for the collection that the cut is acting on. That way other paired collections can access the correct flag for the current cut.
884 else currentObject = objectsToFlag.at(currentObjectIndex);
885 if (currentObjectIndex >= 0 && currentObject == currentCut.inputCollection) continue; // Flags have already been set for the inputCollection object, so should not be set again.
886
887 // single object collections
888 if (currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
889 else if(currentObject == "secondary jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"secondary jets");
890 else if(currentObject == "secondary photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"secondary photons");
891 else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
892 else if(currentObject == "secondary muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"secondary muons");
893 else if(currentObject == "secondary electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"secondary electrons");
894 else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
895 else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
896 else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
897 else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
898 else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
899 else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
900 else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
901 else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
902 else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
903 else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
904 else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
905 else if(currentObject == "trigobjs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,trigobjs.product(),"trigobjs");
906
907
908 // paired object collections
909 else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), "muon-muon pairs");
910 else if(currentObject == "muon-secondary muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), "muon-secondary muon pairs");
911
912 else if(currentObject == "muon-secondary photon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),photons.product(), "muon-secondary photon pairs");
913 else if(currentObject == "muon-secondary jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),jets.product(), "muon-secondary jet pairs");
914 else if(currentObject == "photon-secondary jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),jets.product(), "photon-secondary jet pairs");
915 else if(currentObject == "electron-secondary jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),jets.product(), "electron-secondary jet pairs");
916 else if(currentObject == "electron-secondary electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), "electron-secondary electron pairs");
917
918 else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), "electron-electron pairs");
919 else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), "electron-muon pairs");
920 else if(currentObject == "jet-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),jets.product(), "jet-jet pairs");
921 else if(currentObject == "jet-secondary jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),jets.product(), "jet-secondary jet pairs");
922 else if(currentObject == "electron-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),jets.product(), "electron-jet pairs");
923 else if(currentObject == "electron-photon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),photons.product(), "electron-photon pairs");
924 else if(currentObject == "photon-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),jets.product(), "photon-jet pairs");
925 else if(currentObject == "muon-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),jets.product(), "muon-jet pairs");
926 else if(currentObject == "muon-event pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),events.product(), "muon-event pairs");
927 else if(currentObject == "met-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),jets.product(), "met-jet pairs");
928 else if(currentObject == "track-jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),jets.product(), "track-jet pairs");
929 else if(currentObject == "muon-photon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),photons.product(), "muon-photon pairs");
930 else if(currentObject == "track-event pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),events.product(), "track-event pairs");
931 else if(currentObject == "electron-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),tracks.product(),"electron-track pairs");
932 else if(currentObject == "muon-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),tracks.product(),"muon-track pairs");
933 else if(currentObject == "muon-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),taus.product(),"muon-tau pairs");
934 else if(currentObject == "tau-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),taus.product(),"tau-tau pairs");
935 else if(currentObject == "tau-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus .product(),tracks.product(),"tau-track pairs");
936 else if(currentObject == "electron-trigobj pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),trigobjs.product(),"electron-trigobj pairs");
937 else if(currentObject == "muon-trigobj pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),trigobjs.product(),"muon-trigobj pairs");
938
939 if(currentObject == "stops" && datasetType_ == "signalMC") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,stops.product(),"stops");
940 }
941
942 }//end loop over all cuts
943
944
945
946 //use cumulative flags to apply cuts at event level
947 vector<bool> eventPassedPreviousCuts; //a vector to store cumulative cut descisions after each cut.
948 bool eventPassedAllCuts = true;
949 eventPassedAllCuts = eventPassedAllCuts && triggerDecision; //apply trigger (true if none were specified)
950
951 for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
952
953 //loop over all objects and count how many passed the cumulative selection up to this point
954 cut currentCut = currentChannel.cuts.at(currentCutIndex);
955 int numberPassing = 0;
956 int numberPassingPrev = 0; // number of objects that pass cuts up to the previous one
957
958 for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
959 if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object).first) numberPassing++;
960 }
961 bool cutDecision;
962 if (!currentCut.isVeto) {
963 cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
964 } else {
965 int prevCutIndex = currentCutIndex - 1;
966 if (prevCutIndex<0) {
967 numberPassingPrev = cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size(); // count all objects in collection if cut is the first
968 } else {
969 for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(prevCutIndex).size() ; object++){
970 if (cumulativeFlags.at(currentCut.inputCollection).at(prevCutIndex).at(object).first) {
971 numberPassingPrev++;
972 if (verbose_>1) clog << " object passed previous cut" << endl;
973 }
974 }
975 }
976 int numberFailCut = numberPassingPrev - numberPassing;
977 cutDecision = evaluateComparison(numberFailCut,currentCut.eventComparativeOperator,currentCut.numberRequired);
978 }
979
980 cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
981 eventPassedAllCuts = eventPassedAllCuts && cutDecision;
982 eventPassedPreviousCuts.push_back(eventPassedAllCuts);
983 if (verbose_>1) clog << " Event passed cuts up to cut index " << currentCutIndex << endl;
984
985 }
986 //applying all appropriate scale factors
987 double scaleFactor = masterScaleFactor;
988 muonScaleFactor_ = electronScaleFactor_ = bTagScaleFactor_ = 1.0;
989
990 if(applyLeptonSF_ && datasetType_ != "data"){
991 //only apply SFs if we've cut on this object
992 if(find(objectsToCut.begin(),objectsToCut.end(),"muons") != objectsToCut.end ()){
993 flagPair muonFlags;
994 //get the last valid flags in the flag map
995 for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
996 if (cumulativeFlags.at("muons").at(i).size()){
997 muonFlags = cumulativeFlags.at("muons").at(i);
998 break;
999 }
1000 }
1001 //apply the weight for each of those objects
1002 for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
1003 if(!muonFlags.at(muonIndex).second) continue;
1004 muonScaleFactor_ *= muonSFWeight_->at (muons->at(muonIndex).eta);
1005 }
1006 }
1007 //only apply SFs if we've cut on this object
1008 if(find(objectsToCut.begin(),objectsToCut.end(),"electrons") != objectsToCut.end ()){
1009 flagPair electronFlags;
1010 //get the last valid flags in the flag map
1011 for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
1012 if (cumulativeFlags.at("electrons").at(i).size()){
1013 electronFlags = cumulativeFlags.at("electrons").at(i);
1014 break;
1015 }
1016 }
1017 //apply the weight for each of those objects
1018 for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
1019 if(!electronFlags.at(electronIndex).second) continue;
1020 electronScaleFactor_ *= electronSFWeight_->at (electrons->at(electronIndex).eta, electrons->at(electronIndex).pt);
1021 }
1022 }
1023 }
1024 if(applyBtagSF_ && datasetType_ != "data"){
1025 //only apply SFs if we've cut on this object
1026 if(find(objectsToCut.begin(),objectsToCut.end(),"jets") != objectsToCut.end ()){
1027 flagPair jetFlags;
1028 vector<double> jetSFs;
1029 //get the last valid flags in the flag map
1030 for (int i = cumulativeFlags.at("jets").size() - 1; i >= 0; i--){
1031 if (cumulativeFlags.at("jets").at(i).size()){
1032 jetFlags = cumulativeFlags.at("jets").at(i);
1033 break;
1034 }
1035 }
1036 //apply the weight for each of those objects
1037 for (uint jetIndex = 0; jetIndex != jetFlags.size(); jetIndex++){
1038 if(!jetFlags.at(jetIndex).second) continue;
1039 double jetSFTmp = bTagSFWeight_->sflookup(jets->at(jetIndex).btagCombinedSecVertex, jets->at(jetIndex).pt, jets->at(jetIndex).flavour, jets->at(jetIndex).eta);
1040 jetSFs.push_back(jetSFTmp);
1041 }
1042 bTagScaleFactor_ *= bTagSFWeight_->weight( jetSFs, minBtag_);
1043 }
1044 }
1045 scaleFactor *= muonScaleFactor_;
1046 scaleFactor *= electronScaleFactor_;
1047 scaleFactor *= bTagScaleFactor_;
1048 cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
1049 if (verbose_>1) clog << " Scale factors applied: "
1050 << " muonScaleFactor_ = " << muonScaleFactor_
1051 << ", electronScaleFactor_ = " << electronScaleFactor_
1052 << ", bTagScaleFactor_ = " << bTagScaleFactor_
1053 << ", total scale factor = " << scaleFactor
1054 << endl;
1055
1056
1057 if (printEventInfo_ && eventPassedAllCuts) {
1058 // Write information about event to screen, for testing purposes.
1059 clog << "Event passed all cuts in channel " << currentChannel.name
1060 << ": run:lumi:evt= " << events->at(0).run
1061 << ":" << events->at(0).lumi
1062 << ":" << events->at(0).evt
1063 << endl;
1064 if (findEventsLog) {
1065 (*findEventsLog) << events->at(0).run
1066 << ":" << events->at(0).lumi
1067 << ":" << events->at(0).evt << endl;
1068 }
1069
1070 }
1071
1072
1073 //filling histograms
1074 for(uint currentCut = 0; currentCut != oneDHists_.at(currentChannelIndex).size(); currentCut++){//loop over all the cuts in each channel; if GetPlotsAfterEachCut_ is false, only the last cut will be used.
1075
1076 if (verbose_>2) clog << " Filling histograms for currentcut = " << currentCut << endl;
1077
1078 uint currentDir;
1079 if (!GetPlotsAfterEachCut_) { currentDir = currentChannel.cuts.size() - oneDHists_.at(currentChannelIndex).size(); } //if GetPlotsAfterEachCut_ is false, set currentDir point to the last cut.
1080 else{
1081 currentDir = currentCut;
1082 }
1083
1084
1085 if(eventPassedPreviousCuts.at(currentDir)){
1086
1087 for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
1088 histogram currentHistogram = histograms.at(histogramIndex);
1089
1090 if (cumulativeFlags.count(currentHistogram.inputCollection) == 0) clog << "Error: no flags found for collection: " << currentHistogram.inputCollection << ", will cause a seg fault" << endl;
1091
1092 if (verbose_>1) clog << " Filling histogram " << currentHistogram.name << " for collection " << currentHistogram.inputCollection << endl;
1093
1094 if(currentHistogram.inputVariables.size() == 1){
1095 TH1D* histo;
1096 histo = oneDHists_.at(currentChannelIndex).at(currentCut).at(currentHistogram.name);
1097 if (currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").at(currentDir),scaleFactor);
1098 else if(currentHistogram.inputCollection == "secondary jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("secondary jets").at(currentDir),scaleFactor);
1099 else if(currentHistogram.inputCollection == "secondary photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("secondary photons").at(currentDir),scaleFactor);
1100 else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").at(currentDir),scaleFactor);
1101 else if(currentHistogram.inputCollection == "secondary muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").at(currentDir),scaleFactor);
1102 else if(currentHistogram.inputCollection == "secondary electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("secondary electrons").at(currentDir),scaleFactor);
1103 else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1104 cumulativeFlags.at("muon-muon pairs").at(currentDir),scaleFactor);
1105 else if(currentHistogram.inputCollection == "muon-secondary muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1106 cumulativeFlags.at("muon-secondary muon pairs").at(currentDir),scaleFactor);
1107 else if(currentHistogram.inputCollection == "muon-secondary photon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1108 cumulativeFlags.at("muon-secondary photon pairs").at(currentDir),scaleFactor);
1109 else if(currentHistogram.inputCollection == "muon-secondary jet pairs") fill1DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1110 cumulativeFlags.at("muon-secondary jet pairs").at(currentDir),scaleFactor);
1111 else if(currentHistogram.inputCollection == "photon-secondary jet pairs") fill1DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1112 cumulativeFlags.at("photon-secondary jet pairs").at(currentDir),scaleFactor);
1113 else if(currentHistogram.inputCollection == "electron-secondary jet pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1114 cumulativeFlags.at("electron-secondary jet pairs").at(currentDir),scaleFactor);
1115
1116 else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").at(currentDir),scaleFactor);
1117 else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1118 cumulativeFlags.at("electron-electron pairs").at(currentDir),scaleFactor);
1119 else if(currentHistogram.inputCollection == "jet-jet pairs") fill1DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1120 cumulativeFlags.at("jet-jet pairs").at(currentDir),scaleFactor);
1121 else if(currentHistogram.inputCollection == "jet-secondary jet pairs") fill1DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1122 cumulativeFlags.at("jet-secondary jet pairs").at(currentDir),scaleFactor);
1123
1124 else if(currentHistogram.inputCollection == "electron-secondary electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1125 cumulativeFlags.at("electron-secondary electron pairs").at(currentDir),scaleFactor);
1126 else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(),
1127 cumulativeFlags.at("electron-muon pairs").at(currentDir),scaleFactor);
1128 else if(currentHistogram.inputCollection == "electron-jet pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),jets.product(),
1129 cumulativeFlags.at("electron-jet pairs").at(currentDir),scaleFactor);
1130 else if(currentHistogram.inputCollection == "photon-jet pairs") fill1DHistogram(histo,currentHistogram, photons.product(),jets.product(),
1131 cumulativeFlags.at("photon-jet pairs").at(currentDir),scaleFactor);
1132 else if(currentHistogram.inputCollection == "muon-jet pairs") fill1DHistogram(histo,currentHistogram, muons.product(),jets.product(),
1133 cumulativeFlags.at("muon-jet pairs").at(currentDir),scaleFactor);
1134 else if(currentHistogram.inputCollection == "muon-event pairs") fill1DHistogram(histo,currentHistogram, muons.product(),events.product(),
1135 cumulativeFlags.at("muon-event pairs").at(currentDir),scaleFactor);
1136 else if(currentHistogram.inputCollection == "met-jet pairs") fill1DHistogram(histo,currentHistogram, mets.product(),jets.product(),
1137 cumulativeFlags.at("met-jet pairs").at(currentDir),scaleFactor);
1138 else if(currentHistogram.inputCollection == "track-jet pairs") fill1DHistogram(histo,currentHistogram,tracks.product(),jets.product(),
1139 cumulativeFlags.at("track-jet pairs").at(currentDir),scaleFactor);
1140 else if(currentHistogram.inputCollection == "muon-photon pairs") fill1DHistogram(histo,currentHistogram, muons.product(),photons.product(),
1141 cumulativeFlags.at("muon-photon pairs").at(currentDir),scaleFactor);
1142 else if(currentHistogram.inputCollection == "electron-photon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),photons.product(),
1143 cumulativeFlags.at("electron-photon pairs").at(currentDir),scaleFactor);
1144 else if(currentHistogram.inputCollection == "electron-track pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),tracks.product(),
1145 cumulativeFlags.at("electron-track pairs").at(currentDir),scaleFactor);
1146 else if(currentHistogram.inputCollection == "muon-track pairs") fill1DHistogram(histo,currentHistogram, muons.product(),tracks.product(),
1147 cumulativeFlags.at("muon-track pairs").at(currentDir),scaleFactor);
1148 else if(currentHistogram.inputCollection == "muon-tau pairs") fill1DHistogram(histo,currentHistogram, muons.product(),taus.product(),
1149 cumulativeFlags.at("muon-tau pairs").at(currentDir),scaleFactor);
1150 else if(currentHistogram.inputCollection == "tau-tau pairs") fill1DHistogram(histo,currentHistogram, taus.product(),taus.product(),
1151 cumulativeFlags.at("tau-tau pairs").at(currentDir),scaleFactor);
1152 else if(currentHistogram.inputCollection == "tau-track pairs") fill1DHistogram(histo,currentHistogram, taus.product(),tracks.product(),
1153 cumulativeFlags.at("tau-track pairs").at(currentDir),scaleFactor);
1154 else if(currentHistogram.inputCollection == "electron-trigobj pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),trigobjs.product(),
1155 cumulativeFlags.at("electron-trigobj pairs").at(currentDir),scaleFactor);
1156 else if(currentHistogram.inputCollection == "muon-trigobj pairs") fill1DHistogram(histo,currentHistogram, muons.product(),trigobjs.product(),
1157 cumulativeFlags.at("muon-trigobj pairs").at(currentDir),scaleFactor);
1158
1159 else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").at(currentDir),scaleFactor);
1160 else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").at(currentDir),scaleFactor);
1161 else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").at(currentDir),scaleFactor);
1162 else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").at(currentDir),scaleFactor);
1163 else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").at(currentDir),scaleFactor);
1164 else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").at(currentDir),scaleFactor);
1165 else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").at(currentDir),scaleFactor);
1166 else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").at(currentDir),scaleFactor);
1167 else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").at(currentDir),scaleFactor);
1168 else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").at(currentDir),scaleFactor);
1169 else if(currentHistogram.inputCollection == "trigobjs") fill1DHistogram(histo,currentHistogram,trigobjs.product(),cumulativeFlags.at("trigobjs").at(currentDir),scaleFactor);
1170 else if(currentHistogram.inputCollection == "stops" && datasetType_ == "signalMC") fill1DHistogram(histo,currentHistogram,stops.product(),cumulativeFlags.at("stops").at(currentDir),scaleFactor);
1171 }
1172 else if(currentHistogram.inputVariables.size() == 2){
1173 TH2D* histo;
1174 histo = twoDHists_.at(currentChannelIndex).at(currentCut).at(currentHistogram.name);
1175 if (currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").at(currentDir),scaleFactor);
1176 else if(currentHistogram.inputCollection == "secondary jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("secondary jets").at(currentDir),scaleFactor);
1177 else if(currentHistogram.inputCollection == "secondary photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("secondary photons").at(currentDir),scaleFactor);
1178 else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").at(currentDir),scaleFactor);
1179 else if(currentHistogram.inputCollection == "secondary muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").at(currentDir),scaleFactor);
1180 else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1181 cumulativeFlags.at("muon-muon pairs").at(currentDir),scaleFactor);
1182 else if(currentHistogram.inputCollection == "muon-secondary muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1183 cumulativeFlags.at("muon-secondary muon pairs").at(currentDir),scaleFactor);
1184 else if(currentHistogram.inputCollection == "muon-secondary photon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1185 cumulativeFlags.at("muon-secondary photon pairs").at(currentDir),scaleFactor);
1186 else if(currentHistogram.inputCollection == "muon-secondary jet pairs") fill2DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1187 cumulativeFlags.at("muon-secondary jet pairs").at(currentDir),scaleFactor);
1188 else if(currentHistogram.inputCollection == "electron-secondary jet pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1189 cumulativeFlags.at("electron-secondary jet pairs").at(currentDir),scaleFactor);
1190 else if(currentHistogram.inputCollection == "photon-secondary jet pairs") fill2DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1191 cumulativeFlags.at("photon-secondary jet pairs").at(currentDir),scaleFactor);
1192 else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").at(currentDir),scaleFactor);
1193 else if(currentHistogram.inputCollection == "secondary electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("secondary electrons").at(currentDir),scaleFactor);
1194 else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1195 cumulativeFlags.at("electron-electron pairs").at(currentDir),scaleFactor);
1196 else if(currentHistogram.inputCollection == "jet-jet pairs") fill2DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1197 cumulativeFlags.at("jet-jet pairs").at(currentDir),scaleFactor);
1198 else if(currentHistogram.inputCollection == "electron-secondary electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1199 cumulativeFlags.at("electron-secondary electron pairs").at(currentDir),scaleFactor);
1200 else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(),
1201 cumulativeFlags.at("electron-muon pairs").at(currentDir),scaleFactor);
1202 else if(currentHistogram.inputCollection == "electron-jet pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1203 cumulativeFlags.at("electron-jet pairs").at(currentDir),scaleFactor);
1204 else if(currentHistogram.inputCollection == "electron-photon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),photons.product(),
1205 cumulativeFlags.at("electron-photon pairs").at(currentDir),scaleFactor);
1206 else if(currentHistogram.inputCollection == "muon-jet pairs") fill2DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1207 cumulativeFlags.at("muon-jet pairs").at(currentDir),scaleFactor);
1208 else if(currentHistogram.inputCollection == "muon-event pairs") fill2DHistogram(histo,currentHistogram,muons.product(),events.product(),
1209 cumulativeFlags.at("muon-event pairs").at(currentDir),scaleFactor);
1210 else if(currentHistogram.inputCollection == "met-jet pairs") fill2DHistogram(histo,currentHistogram,mets.product(),jets.product(),
1211 cumulativeFlags.at("met-jet pairs").at(currentDir),scaleFactor);
1212 else if(currentHistogram.inputCollection == "track-jet pairs") fill2DHistogram(histo,currentHistogram,tracks.product(),jets.product(),
1213 cumulativeFlags.at("track-jet pairs").at(currentDir),scaleFactor);
1214 else if(currentHistogram.inputCollection == "photon-jet pairs") fill2DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1215 cumulativeFlags.at("photon-jet pairs").at(currentDir),scaleFactor);
1216 else if(currentHistogram.inputCollection == "muon-photon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1217 cumulativeFlags.at("muon-photon pairs").at(currentDir),scaleFactor);
1218 else if(currentHistogram.inputCollection == "electron-track pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),tracks.product(),
1219 cumulativeFlags.at("electron-track pairs").at(currentDir),scaleFactor);
1220 else if(currentHistogram.inputCollection == "muon-track pairs") fill2DHistogram(histo,currentHistogram,muons.product(),tracks.product(),
1221 cumulativeFlags.at("muon-track pairs").at(currentDir),scaleFactor);
1222 else if(currentHistogram.inputCollection == "muon-tau pairs") fill2DHistogram(histo,currentHistogram,muons.product(),taus.product(),
1223 cumulativeFlags.at("muon-tau pairs").at(currentDir),scaleFactor);
1224 else if(currentHistogram.inputCollection == "tau-tau pairs") fill2DHistogram(histo,currentHistogram,taus.product(),taus.product(),
1225 cumulativeFlags.at("tau-tau pairs").at(currentDir),scaleFactor);
1226 else if(currentHistogram.inputCollection == "tau-track pairs") fill2DHistogram(histo,currentHistogram,taus.product(),tracks.product(),
1227 cumulativeFlags.at("tau-track pairs").at(currentDir),scaleFactor);
1228 else if(currentHistogram.inputCollection == "electron-trigobj pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),trigobjs.product(),
1229 cumulativeFlags.at("electron-trigobj pairs").at(currentDir),scaleFactor);
1230 else if(currentHistogram.inputCollection == "muon-trigobj pairs") fill2DHistogram(histo,currentHistogram,muons.product(),trigobjs.product(),
1231 cumulativeFlags.at("muon-trigobj pairs").at(currentDir),scaleFactor);
1232 else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").at(currentDir),scaleFactor);
1233 else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").at(currentDir),scaleFactor);
1234 else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").at(currentDir),scaleFactor);
1235 else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").at(currentDir),scaleFactor);
1236 else if(currentHistogram.inputCollection == "track-event pairs") fill2DHistogram(histo,currentHistogram,tracks.product(),events.product(),
1237 cumulativeFlags.at("track-event pairs").at(currentDir),scaleFactor);
1238 else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").at(currentDir),scaleFactor);
1239 else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").at(currentDir),scaleFactor);
1240 else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").at(currentDir),scaleFactor);
1241 else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").at(currentDir),scaleFactor);
1242 else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").at(currentDir),scaleFactor);
1243 else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").at(currentDir),scaleFactor);
1244 else if(currentHistogram.inputCollection == "trigobjs") fill2DHistogram(histo,currentHistogram,trigobjs.product(),cumulativeFlags.at("trigobjs").at(currentDir),scaleFactor);
1245 else if(currentHistogram.inputCollection == "stops" && datasetType_ == "signalMC") fill2DHistogram(histo,currentHistogram,stops.product(),cumulativeFlags.at("stops").at(currentDir),scaleFactor);
1246 }
1247
1248 }
1249
1250
1251 //fills histograms with the sizes of collections
1252 for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
1253
1254 string currentObject = objectsToPlot.at(currentObjectIndex);
1255 string objectToPlot = "";
1256
1257 if (verbose_) clog << "Filling histogram of number of selected objects in collection: " << currentObject << endl;
1258
1259 // Name of objectToPlot here must match the name specified in OSUAnalysis::OSUAnalysis().
1260 if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
1261 else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
1262 else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
1263 else if(currentObject == "electron-photon pairs") objectToPlot = "electronPhotonPairs";
1264 else if(currentObject == "electron-jet pairs") objectToPlot = "electronJetPairs";
1265 else if(currentObject == "muon-jet pairs") objectToPlot = "muonJetPairs";
1266 else if(currentObject == "muon-event pairs") objectToPlot = "muonEventPairs";
1267 else if(currentObject == "muon-photon pairs") objectToPlot = "muonPhotonPairs";
1268 else if(currentObject == "photon-jet pairs") objectToPlot = "photonJetPairs";
1269 else if(currentObject == "met-jet pairs") objectToPlot = "metJetPairs";
1270 else if(currentObject == "track-jet pairs") objectToPlot = "trackJetPairs";
1271 else if(currentObject == "jet-jet pairs") objectToPlot = "dijetPairs";
1272 else if(currentObject == "jet-secondary jet pairs") objectToPlot = "jetSecondaryJetPairs";
1273 else if(currentObject == "secondary jets") objectToPlot = "secondaryJets";
1274 else if(currentObject == "secondary photons") objectToPlot = "secondaryPhotons";
1275 else if(currentObject == "electron-track pairs") objectToPlot = "electronTrackPairs";
1276 else if(currentObject == "muon-track pairs") objectToPlot = "muonTrackPairs";
1277 else if(currentObject == "muon-tau pairs") objectToPlot = "muonTauPairs";
1278 else if(currentObject == "tau-tau pairs") objectToPlot = "ditauPairs";
1279 else if(currentObject == "tau-track pairs") objectToPlot = "tauTrackPairs";
1280 else if(currentObject == "track-event pairs") objectToPlot = "trackEventPairs";
1281 else if(currentObject == "muon-secondary muon pairs") objectToPlot = "muonSecondaryMuonPairs";
1282 else if(currentObject == "secondary muons") objectToPlot = "secondaryMuons";
1283 else if(currentObject == "muon-secondary jet pairs") objectToPlot = "muonSecondaryJetPairs";
1284 else if(currentObject == "muon-secondary photon pairs") objectToPlot = "muonSecondaryJetPairs";
1285 else if(currentObject == "electron-secondary jet pairs") objectToPlot = "electronSecondaryJetPairs";
1286 else if(currentObject == "photon-secondary jet pairs") objectToPlot = "photonSecondaryJetPairs";
1287 else if(currentObject == "electron-secondary electron pairs") objectToPlot = "electronSecondaryElectronPairs";
1288 else if(currentObject == "secondary electrons") objectToPlot = "secondaryElectrons";
1289 else if(currentObject == "electron-trigobj pairs") objectToPlot = "electronTrigobjPairs";
1290 else if(currentObject == "muon-trigobj pairs") objectToPlot = "muonTrigobjPairs";
1291 else objectToPlot = currentObject;
1292
1293 string tempCurrentObject = objectToPlot;
1294 tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
1295 string histoName = "num" + tempCurrentObject;
1296
1297
1298 if(find(objectsToPlot.begin(), objectsToPlot.end(), currentObject) != objectsToPlot.end()) {
1299 flagPair lastCutFlags = cumulativeFlags.at(currentObject).at(currentDir);
1300 int numToPlot = 0;
1301 for (uint iObj = 0; iObj != lastCutFlags.size(); iObj++){ // loop over all the objects
1302 if(lastCutFlags.at(iObj).second) {
1303 numToPlot++;
1304 if (verbose_>3) clog << " Found object " << iObj << " in collection " << currentObject << " to plot." << endl;
1305 }
1306 }
1307
1308 if(objectToPlot == "primaryvertexs"){
1309 oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
1310 oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
1311 }
1312 else {
1313 if (printEventInfo_) clog << "Number of selected " << objectToPlot << " to plot: " << numToPlot << endl;
1314 oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName)->Fill(numToPlot,scaleFactor);
1315 }
1316 }
1317
1318 } // end for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++)
1319
1320
1321 } // end if(eventPassedPreviousCuts.at(currentDir))
1322 } // end loop over cuts
1323
1324 if (!useEDMFormat_ && eventPassedAllCuts){
1325 // Assign BNTree variables
1326 for (uint iBranch = 0; iBranch != treeBranches_.size(); iBranch++) {
1327 BranchSpecs brSpecs = treeBranches_.at(iBranch);
1328 string coll = brSpecs.inputCollection;
1329 if (cumulativeFlags.count(coll) == 0) clog << "Error: no flags found for collection: " << coll << ", will cause a seg fault" << endl;
1330
1331 if (coll == "jets") assignTreeBranch(brSpecs,jets.product(), cumulativeFlags.at(coll).back());
1332 else if(coll == "secondary jets") assignTreeBranch(brSpecs,jets.product(), cumulativeFlags.at(coll).back());
1333 else if(coll == "secondary photons") assignTreeBranch(brSpecs,photons.product(), cumulativeFlags.at(coll).back());
1334 else if(coll == "muons") assignTreeBranch(brSpecs,muons.product(), cumulativeFlags.at(coll).back());
1335 else if(coll == "secondary muons") assignTreeBranch(brSpecs,muons.product(), cumulativeFlags.at(coll).back());
1336 else if(coll == "electrons") assignTreeBranch(brSpecs,electrons.product(), cumulativeFlags.at(coll).back());
1337 else if(coll == "secondary electrons") assignTreeBranch(brSpecs,electrons.product(), cumulativeFlags.at(coll).back());
1338 else if(coll == "events") assignTreeBranch(brSpecs,events.product(), cumulativeFlags.at(coll).back());
1339 else if(coll == "taus") assignTreeBranch(brSpecs,taus.product(), cumulativeFlags.at(coll).back());
1340 else if(coll == "mets") assignTreeBranch(brSpecs,mets.product(), cumulativeFlags.at(coll).back());
1341 else if(coll == "tracks") assignTreeBranch(brSpecs,tracks.product(), cumulativeFlags.at(coll).back());
1342 else if(coll == "genjets") assignTreeBranch(brSpecs,genjets.product(), cumulativeFlags.at(coll).back());
1343 else if(coll == "mcparticles") assignTreeBranch(brSpecs,mcparticles.product(), cumulativeFlags.at(coll).back());
1344 else if(coll == "primaryvertexs") assignTreeBranch(brSpecs,primaryvertexs.product(),cumulativeFlags.at(coll).back());
1345 else if(coll == "bxlumis") assignTreeBranch(brSpecs,bxlumis.product(), cumulativeFlags.at(coll).back());
1346 else if(coll == "photons") assignTreeBranch(brSpecs,photons.product(), cumulativeFlags.at(coll).back());
1347 else if(coll == "superclusters") assignTreeBranch(brSpecs,superclusters.product(), cumulativeFlags.at(coll).back());
1348 else if(coll == "trigobjs") assignTreeBranch(brSpecs,trigobjs.product(), cumulativeFlags.at(coll).back());
1349 else if(coll == "stops"
1350 && datasetType_ == "signalMC") assignTreeBranch(brSpecs,stops.product(), cumulativeFlags.at(coll).back());
1351 } // end loop over branches
1352 // set the evtLong, runInt, and lumiInt variables separately, since they are not type float
1353 BNTreeBranchVals_evtLong_ = events->at(0).evt;
1354 BNTreeBranchVals_runInt_ = events->at(0).run;
1355 BNTreeBranchVals_lumiInt_ = events->at(0).lumi;
1356
1357 if (!BNTrees_.at(currentChannelIndex)) { clog << "ERROR: Undefined BNTree. Will likely seg fault." << endl; }
1358 BNTrees_.at(currentChannelIndex)->Fill(); // only fill if event has passed cuts
1359 }
1360
1361 (*channelDecisions)[currentChannel.name] = eventPassedAllCuts;
1362
1363 } // end loop over channel
1364
1365 masterCutFlow_->fillCutFlow(masterScaleFactor);
1366
1367 event.put (channelDecisions, "channelDecisions");
1368
1369 isFirstEvent_ = false;
1370
1371 if (verbose_) clog << "Finished OSUAnalysis::produce." << endl;
1372
1373 } // end void OSUAnalysis::produce (const edm::Event &event, const edm::EventSetup &setup)
1374
1375
1376
1377 bool
1378 OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
1379
1380
1381 if(comparison == ">") return testValue > cutValue;
1382 else if(comparison == ">=") return testValue >= cutValue;
1383 else if(comparison == "<") return testValue < cutValue;
1384 else if(comparison == "<=") return testValue <= cutValue;
1385 else if(comparison == "==") return testValue == cutValue;
1386 else if(comparison == "=") return testValue == cutValue;
1387 else if(comparison == "!=") return testValue != cutValue;
1388 else {clog << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
1389
1390 }
1391
1392 bool
1393 OSUAnalysis::evaluateComparison (string testValue, string comparison, string cutValue){
1394
1395
1396 if(comparison == ">") return testValue > cutValue;
1397 else if(comparison == ">=") return testValue >= cutValue;
1398 else if(comparison == "<") return testValue < cutValue;
1399 else if(comparison == "<=") return testValue <= cutValue;
1400 else if(comparison == "==") return testValue == cutValue;
1401 else if(comparison == "=") return testValue == cutValue;
1402 else if(comparison == "!=") return testValue != cutValue;
1403 else {clog << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
1404
1405 }
1406
1407 bool
1408 OSUAnalysis::evaluateTriggers (vector<string> triggersToTest, vector<string> triggersToVeto, const BNtriggerCollection* triggerCollection){
1409
1410 //initialize to false until a chosen trigger is passed
1411 bool triggerDecision = false;
1412
1413 if (printAllTriggers_) clog << "Printing list of all available triggers (which this event may or may not pass):" << endl;
1414 //loop over all triggers defined in the event
1415 for (BNtriggerCollection::const_iterator trigger = triggerCollection->begin (); trigger != triggerCollection->end (); trigger++){
1416
1417 if (printAllTriggers_) clog << " " << trigger->name << endl;
1418
1419 //we're only interested in triggers that actually passed
1420 if(trigger->pass != 1) continue;
1421
1422 //if the event passes one of the veto triggers, exit and return false
1423 for(uint triggerName = 0; triggerName != triggersToVeto.size(); triggerName++){
1424 if(trigger->name.find(triggersToVeto.at(triggerName))!=string::npos) return false;
1425 }
1426 //if the event passes one of the chosen triggers, set triggerDecision to true
1427 for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++){
1428 if(trigger->name.find(triggersToTest.at(triggerName))!=string::npos) triggerDecision = true;
1429 }
1430 }
1431
1432 printAllTriggers_ = false; // only print triggers once, not every event
1433
1434 //if none of the veto triggers fired:
1435 //return the OR of all the chosen triggers
1436 if (triggersToTest.size() != 0) return triggerDecision;
1437 //or if no triggers were defined return true
1438 else return true;
1439 }
1440
1441
1442 vector<string>
1443 OSUAnalysis::splitString (string inputString){
1444
1445 stringstream stringStream(inputString);
1446 istream_iterator<string> begin(stringStream);
1447 istream_iterator<string> end;
1448 vector<string> stringVector(begin, end);
1449 return stringVector;
1450
1451 }
1452
1453
1454 void OSUAnalysis::getTwoObjs(string tempInputCollection, string& obj1, string& obj2) {
1455 // Set two object strings from the tempInputCollection string,
1456 // For example, if tempInputCollection is "electron-muon pairs",
1457 // then obj1 = "electrons" and obj2 = "muons".
1458 // Note that the objects have an "s" appended.
1459
1460 int dashIndex = tempInputCollection.find("-");
1461 int spaceIndex = tempInputCollection.find_last_of(" ");
1462 int secondWordLength = spaceIndex - dashIndex;
1463 obj1 = tempInputCollection.substr(0,dashIndex) + "s";
1464 obj2 = tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s";
1465
1466 }
1467
1468
1469 string OSUAnalysis::getObjToGet(string obj) {
1470 // Return the string corresponding to the object to get for the given obj string.
1471 // Right now this only handles the case in which obj contains "secondary",
1472 // e.g, "secondary muons".
1473 // Note that "s" is NOT appended.
1474
1475 if (obj.find("secondary")==string::npos) return obj; // "secondary" is not found
1476 int firstSpaceIndex = obj.find_first_of(" ");
1477 return obj.substr(firstSpaceIndex+1,obj.length()-1);
1478
1479 }
1480
1481
1482 //!jet valueLookup
1483 double
1484 OSUAnalysis::valueLookup (const BNjet* object, string variable, string function, string &stringValue){
1485
1486 double value = 0.0;
1487 if(variable == "energy") value = object->energy;
1488 else if(variable == "et") value = object->et;
1489 else if(variable == "pt") value = object->pt;
1490 else if(variable == "px") value = object->px;
1491 else if(variable == "py") value = object->py;
1492 else if(variable == "pz") value = object->pz;
1493 else if(variable == "phi") value = object->phi;
1494 else if(variable == "eta") value = object->eta;
1495 else if(variable == "theta") value = object->theta;
1496 else if(variable == "Upt") value = object->Upt;
1497 else if(variable == "Uenergy") value = object->Uenergy;
1498 else if(variable == "L2pt") value = object->L2pt;
1499 else if(variable == "L2L3pt") value = object->L2L3pt;
1500 else if(variable == "L2L3respt") value = object->L2L3respt;
1501 else if(variable == "respt") value = object->respt;
1502 else if(variable == "EMfrac") value = object->EMfrac;
1503 else if(variable == "Hadfrac") value = object->Hadfrac;
1504 else if(variable == "charge") value = object->charge;
1505 else if(variable == "mass") value = object->mass;
1506 else if(variable == "area") value = object->area;
1507 else if(variable == "fHPD") value = object->fHPD;
1508 else if(variable == "approximatefHPD") value = object->approximatefHPD;
1509 else if(variable == "genPartonET") value = object->genPartonET;
1510 else if(variable == "genPartonPT") value = object->genPartonPT;
1511 else if(variable == "genPartonEta") value = object->genPartonEta;
1512 else if(variable == "genPartonPhi") value = object->genPartonPhi;
1513 else if(variable == "genJetET") value = object->genJetET;
1514 else if(variable == "genJetPT") value = object->genJetPT;
1515 else if(variable == "genJetEta") value = object->genJetEta;
1516 else if(variable == "genJetPhi") value = object->genJetPhi;
1517 else if(variable == "btagTChighPur") value = object->btagTChighPur;
1518 else if(variable == "btagTChighEff") value = object->btagTChighEff;
1519 else if(variable == "btagJetProb") value = object->btagJetProb;
1520 else if(variable == "btagJetBProb") value = object->btagJetBProb;
1521 else if(variable == "btagSoftEle") value = object->btagSoftEle;
1522 else if(variable == "btagSoftMuon") value = object->btagSoftMuon;
1523 else if(variable == "btagSoftMuonNoIP") value = object->btagSoftMuonNoIP;
1524 else if(variable == "btagSecVertex") value = object->btagSecVertex;
1525 else if(variable == "btagSecVertexHighEff") value = object->btagSecVertexHighEff;
1526 else if(variable == "btagSecVertexHighPur") value = object->btagSecVertexHighPur;
1527 else if(variable == "btagCombinedSecVertex") value = object->btagCombinedSecVertex;
1528 else if(variable == "btagCombinedSecVertexMVA") value = object->btagCombinedSecVertexMVA;
1529 else if(variable == "btagSoftMuonByPt") value = object->btagSoftMuonByPt;
1530 else if(variable == "btagSoftMuonByIP3") value = object->btagSoftMuonByIP3;
1531 else if(variable == "btagSoftElectronByPt") value = object->btagSoftElectronByPt;
1532 else if(variable == "btagSoftElectronByIP3") value = object->btagSoftElectronByIP3;
1533 else if(variable == "n90Hits") value = object->n90Hits;
1534 else if(variable == "hitsInN90") value = object->hitsInN90;
1535 else if(variable == "chargedHadronEnergyFraction") value = object->chargedHadronEnergyFraction;
1536 else if(variable == "neutralHadronEnergyFraction") value = object->neutralHadronEnergyFraction;
1537 else if(variable == "chargedEmEnergyFraction") value = object->chargedEmEnergyFraction;
1538 else if(variable == "neutralEmEnergyFraction") value = object->neutralEmEnergyFraction;
1539 else if(variable == "fLong") value = object->fLong;
1540 else if(variable == "fShort") value = object->fShort;
1541 else if(variable == "etaetaMoment") value = object->etaetaMoment;
1542 else if(variable == "phiphiMoment") value = object->phiphiMoment;
1543 else if(variable == "JESunc") value = object->JESunc;
1544 else if(variable == "JECuncUp") value = object->JECuncUp;
1545 else if(variable == "JECuncDown") value = object->JECuncDown;
1546 else if(variable == "puJetMVA_full") value = object->puJetMVA_full;
1547 else if(variable == "puJetMVA_simple") value = object->puJetMVA_simple;
1548 else if(variable == "puJetMVA_cutbased") value = object->puJetMVA_cutbased;
1549 else if(variable == "dZ") value = object->dZ;
1550 else if(variable == "dR2Mean") value = object->dR2Mean;
1551 else if(variable == "dRMean") value = object->dRMean;
1552 else if(variable == "frac01") value = object->frac01;
1553 else if(variable == "frac02") value = object->frac02;
1554 else if(variable == "frac03") value = object->frac03;
1555 else if(variable == "frac04") value = object->frac04;
1556 else if(variable == "frac05") value = object->frac05;
1557 else if(variable == "frac06") value = object->frac06;
1558 else if(variable == "frac07") value = object->frac07;
1559 else if(variable == "beta") value = object->beta;
1560 else if(variable == "betaStar") value = object->betaStar;
1561 else if(variable == "betaClassic") value = object->betaClassic;
1562 else if(variable == "betaStarClassic") value = object->betaStarClassic;
1563 else if(variable == "ptD") value = object->ptD;
1564 else if(variable == "nvtx") value = object->nvtx;
1565 else if(variable == "d0") value = object->d0;
1566 else if(variable == "leadCandPt") value = object->leadCandPt;
1567 else if(variable == "leadCandVx") value = object->leadCandVx;
1568 else if(variable == "leadCandVy") value = object->leadCandVy;
1569 else if(variable == "leadCandVz") value = object->leadCandVz;
1570 else if(variable == "leadCandDistFromPV") value = object->leadCandDistFromPV;
1571 else if(variable == "flavour") value = object->flavour;
1572 else if(variable == "Nconst") value = object->Nconst;
1573 else if(variable == "jetIDMinimal") value = object->jetIDMinimal;
1574 else if(variable == "jetIDLooseAOD") value = object->jetIDLooseAOD;
1575 else if(variable == "jetIDLoose") value = object->jetIDLoose;
1576 else if(variable == "jetIDTight") value = object->jetIDTight;
1577 else if(variable == "genPartonId") value = object->genPartonId;
1578 else if(variable == "genPartonMotherId") value = object->genPartonMotherId;
1579 else if(variable == "genPartonMother0Id") value = object->genPartonMother0Id;
1580 else if(variable == "genPartonMother1Id") value = object->genPartonMother1Id;
1581 else if(variable == "genPartonGrandMotherId") value = object->genPartonGrandMotherId;
1582 else if(variable == "genPartonGrandMother00Id") value = object->genPartonGrandMother00Id;
1583 else if(variable == "genPartonGrandMother01Id") value = object->genPartonGrandMother01Id;
1584 else if(variable == "genPartonGrandMother10Id") value = object->genPartonGrandMother10Id;
1585 else if(variable == "genPartonGrandMother11Id") value = object->genPartonGrandMother11Id;
1586 else if(variable == "chargedMultiplicity") value = object->chargedMultiplicity;
1587 else if(variable == "neutralMultiplicity") value = object->neutralMultiplicity;
1588 else if(variable == "nconstituents") value = object->nconstituents;
1589 else if(variable == "nHit") value = object->nHit;
1590 else if(variable == "puJetId_full") value = object->puJetId_full;
1591 else if(variable == "puJetId_simple") value = object->puJetId_simple;
1592 else if(variable == "puJetId_cutbased") value = object->puJetId_cutbased;
1593 else if(variable == "puJetId_tight_full") value = object->puJetId_tight_full;
1594 else if(variable == "puJetId_tight_simple") value = object->puJetId_tight_simple;
1595 else if(variable == "puJetId_tight_cutbased") value = object->puJetId_tight_cutbased;
1596 else if(variable == "puJetId_medium_full") value = object->puJetId_medium_full;
1597 else if(variable == "puJetId_medium_simple") value = object->puJetId_medium_simple;
1598 else if(variable == "puJetId_medium_cutbased") value = object->puJetId_medium_cutbased;
1599 else if(variable == "puJetId_loose_full") value = object->puJetId_loose_full;
1600 else if(variable == "puJetId_loose_simple") value = object->puJetId_loose_simple;
1601 else if(variable == "puJetId_loose_cutbased") value = object->puJetId_loose_cutbased;
1602
1603 //user defined variable
1604 else if(variable == "disappTrkLeadingJetID") {
1605 value = object->pt > 110
1606 && fabs(object->eta) < 2.4
1607 && object->chargedHadronEnergyFraction > 0.2
1608 && object->neutralHadronEnergyFraction < 0.7
1609 && object->chargedEmEnergyFraction < 0.5
1610 && object->neutralEmEnergyFraction < 0.7;
1611 }
1612
1613 else if(variable == "disappTrkSubLeadingJetID") {
1614 value = object->pt > 30
1615 && fabs(object->eta) < 4.5
1616 && object->neutralHadronEnergyFraction < 0.7
1617 && object->chargedEmEnergyFraction < 0.5;
1618 }
1619
1620 else if(variable == "dPhiMet") {
1621 if (const BNmet *met = chosenMET ()) {
1622 value = deltaPhi (object->phi, met->phi);
1623 } else value = -999;
1624 }
1625
1626
1627 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1628
1629 value = applyFunction(function, value);
1630
1631 return value;
1632 } // end jet valueLookup
1633
1634
1635 //!muon valueLookup
1636 double
1637 OSUAnalysis::valueLookup (const BNmuon* object, string variable, string function, string &stringValue){
1638
1639 double value = 0.0;
1640 if(variable == "energy") value = object->energy;
1641 else if(variable == "et") value = object->et;
1642 else if(variable == "pt") value = object->pt;
1643 else if(variable == "px") value = object->px;
1644 else if(variable == "py") value = object->py;
1645 else if(variable == "pz") value = object->pz;
1646 else if(variable == "phi") value = object->phi;
1647 else if(variable == "eta") value = object->eta;
1648 else if(variable == "theta") value = object->theta;
1649 else if(variable == "trackIso") value = object->trackIso;
1650 else if(variable == "ecalIso") value = object->ecalIso;
1651 else if(variable == "hcalIso") value = object->hcalIso;
1652 else if(variable == "caloIso") value = object->caloIso;
1653 else if(variable == "trackIsDR03") value = object->trackIsoDR03;
1654 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
1655 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
1656 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
1657 else if(variable == "trackVetoIsoDR03") value = object->trackVetoIsoDR03;
1658 else if(variable == "ecalVetoIsoDR03") value = object->ecalVetoIsoDR03;
1659 else if(variable == "hcalVetoIsoDR03") value = object->hcalVetoIsoDR03;
1660 else if(variable == "caloVetoIsoDR03") value = object->caloVetoIsoDR03;
1661 else if(variable == "trackIsoDR05") value = object->trackIsoDR05;
1662 else if(variable == "ecalIsoDR05") value = object->ecalIsoDR05;
1663 else if(variable == "hcalIsoDR05") value = object->hcalIsoDR05;
1664 else if(variable == "caloIsoDR05") value = object->caloIsoDR05;
1665 else if(variable == "trackVetoIsoDR05") value = object->trackVetoIsoDR05;
1666 else if(variable == "ecalVetoIsoDR05") value = object->ecalVetoIsoDR05;
1667 else if(variable == "hcalVetoIsoDR05") value = object->hcalVetoIsoDR05;
1668 else if(variable == "caloVetoIsoDR05") value = object->caloVetoIsoDR05;
1669 else if(variable == "hcalE") value = object->hcalE;
1670 else if(variable == "ecalE") value = object->ecalE;
1671 else if(variable == "genET") value = object->genET;
1672 else if(variable == "genPT") value = object->genPT;
1673 else if(variable == "genPhi") value = object->genPhi;
1674 else if(variable == "genEta") value = object->genEta;
1675 else if(variable == "genMotherET") value = object->genMotherET;
1676 else if(variable == "genMotherPT") value = object->genMotherPT;
1677 else if(variable == "genMotherPhi") value = object->genMotherPhi;
1678 else if(variable == "genMotherEta") value = object->genMotherEta;
1679 else if(variable == "vx") value = object->vx;
1680 else if(variable == "vy") value = object->vy;
1681 else if(variable == "vz") value = object->vz;
1682 else if(variable == "tkNormChi2") value = object->tkNormChi2;
1683 else if(variable == "tkPT") value = object->tkPT;
1684 else if(variable == "tkEta") value = object->tkEta;
1685 else if(variable == "tkPhi") value = object->tkPhi;
1686 else if(variable == "tkDZ") value = object->tkDZ;
1687 else if(variable == "tkD0") value = object->tkD0;
1688 else if(variable == "tkD0bs") value = object->tkD0bs;
1689 else if(variable == "tkD0err") value = object->tkD0err;
1690 else if(variable == "samNormChi2") value = object->samNormChi2;
1691 else if(variable == "samPT") value = object->samPT;
1692 else if(variable == "samEta") value = object->samEta;
1693 else if(variable == "samPhi") value = object->samPhi;
1694 else if(variable == "samDZ") value = object->samDZ;
1695 else if(variable == "samD0") value = object->samD0;
1696 else if(variable == "samD0bs") value = object->samD0bs;
1697 else if(variable == "samD0err") value = object->samD0err;
1698 else if(variable == "comNormChi2") value = object->comNormChi2;
1699 else if(variable == "comPT") value = object->comPT;
1700 else if(variable == "comEta") value = object->comEta;
1701 else if(variable == "comPhi") value = object->comPhi;
1702 else if(variable == "comDZ") value = object->comDZ;
1703 else if(variable == "comD0") value = object->comD0;
1704 else if(variable == "comD0bs") value = object->comD0bs;
1705 else if(variable == "comD0err") value = object->comD0err;
1706 else if(variable == "isolationR03emVetoEt") value = object->isolationR03emVetoEt;
1707 else if(variable == "isolationR03hadVetoEt") value = object->isolationR03hadVetoEt;
1708 else if(variable == "normalizedChi2") value = object->normalizedChi2;
1709 else if(variable == "dVzPVz") value = object->dVzPVz;
1710 else if(variable == "dB") value = object->dB;
1711 else if(variable == "ptErr") value = object->ptErr;
1712 else if(variable == "innerTrackNormChi2") value = object->innerTrackNormChi2;
1713 else if(variable == "correctedD0") value = object->correctedD0;
1714 else if(variable == "correctedD0Vertex") value = object->correctedD0Vertex;
1715 else if(variable == "correctedDZ") value = object->correctedDZ;
1716 else if(variable == "particleIso") value = object->particleIso;
1717 else if(variable == "chargedHadronIso") value = object->chargedHadronIso;
1718 else if(variable == "neutralHadronIso") value = object->neutralHadronIso;
1719 else if(variable == "photonIso") value = object->photonIso;
1720 else if(variable == "puChargedHadronIso") value = object->puChargedHadronIso;
1721 else if(variable == "chargedHadronIsoDR03") value = object->chargedHadronIsoDR03;
1722 else if(variable == "neutralHadronIsoDR03") value = object->neutralHadronIsoDR03;
1723 else if(variable == "photonIsoDR03") value = object->photonIsoDR03;
1724 else if(variable == "puChargedHadronIsoDR03") value = object->puChargedHadronIsoDR03;
1725 else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
1726 else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1727 else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1728 else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1729 else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
1730 else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1731 else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1732 else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1733 else if(variable == "rhoPrime") value = object->rhoPrime;
1734 else if(variable == "AEffDr03") value = object->AEffDr03;
1735 else if(variable == "AEffDr04") value = object->AEffDr04;
1736 else if(variable == "pfIsoR03SumChargedHadronPt") value = object->pfIsoR03SumChargedHadronPt;
1737 else if(variable == "pfIsoR03SumNeutralHadronEt") value = object->pfIsoR03SumNeutralHadronEt;
1738 else if(variable == "pfIsoR03SumPhotonEt") value = object->pfIsoR03SumPhotonEt;
1739 else if(variable == "pfIsoR03SumPUPt") value = object->pfIsoR03SumPUPt;
1740 else if(variable == "relpfIsoR04SumExceptChargedHad") value = (object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)/ object->pt;
1741 else if(variable == "relpfIsoR04SumChargedHadronPt") value = object->pfIsoR04SumChargedHadronPt/object->pt;
1742 else if(variable == "relpfIsoR04SumNeutralHadronEt") value = object->pfIsoR04SumNeutralHadronEt/object->pt;
1743 else if(variable == "relpfIsoR04SumPhotonEt") value = object->pfIsoR04SumPhotonEt/object->pt;
1744 else if(variable == "relpfIsoR04SumPUPt") value = object->pfIsoR04SumPUPt/object->pt;
1745 else if(variable == "pfIsoR04SumChargedHadronPt") value = object->pfIsoR04SumChargedHadronPt;
1746 else if(variable == "pfIsoR04SumNeutralHadronEt") value = object->pfIsoR04SumNeutralHadronEt;
1747 else if(variable == "pfIsoR04SumPhotonEt") value = object->pfIsoR04SumPhotonEt;
1748 else if(variable == "pfIsoR04SumPUPt") value = object->pfIsoR04SumPUPt;
1749 else if(variable == "IP") value = object->IP;
1750 else if(variable == "IPError") value = object->IPError;
1751 else if(variable == "timeAtIpInOut") value = object->timeAtIpInOut;
1752 else if(variable == "timeAtIpInOutErr") value = object->timeAtIpInOutErr;
1753 else if(variable == "timeAtIpOutIn") value = object->timeAtIpOutIn;
1754 else if(variable == "timeAtIpOutInErr") value = object->timeAtIpOutInErr;
1755 else if(variable == "ecal_time") value = object->ecal_time;
1756 else if(variable == "hcal_time") value = object->hcal_time;
1757 else if(variable == "ecal_timeError") value = object->ecal_timeError;
1758 else if(variable == "hcal_timeError") value = object->hcal_timeError;
1759 else if(variable == "energy_ecal") value = object->energy_ecal;
1760 else if(variable == "energy_hcal") value = object->energy_hcal;
1761 else if(variable == "e3x3_ecal") value = object->e3x3_ecal;
1762 else if(variable == "e3x3_hcal") value = object->e3x3_hcal;
1763 else if(variable == "energyMax_ecal") value = object->energyMax_ecal;
1764 else if(variable == "energyMax_hcal") value = object->energyMax_hcal;
1765 else if(variable == "charge") value = object->charge;
1766 else if(variable == "IDGMPTight") value = object->IDGMPTight;
1767 else if(variable == "tkNumValidHits") value = object->tkNumValidHits;
1768 else if(variable == "tkCharge") value = object->tkCharge;
1769 else if(variable == "samNumValidHits") value = object->samNumValidHits;
1770 else if(variable == "samCharge") value = object->samCharge;
1771 else if(variable == "comNumValidHits") value = object->comNumValidHits;
1772 else if(variable == "comCharge") value = object->comCharge;
1773 else if(variable == "genId") value = object->genId;
1774 else if(variable == "genCharge") value = object->genCharge;
1775 else if(variable == "genNumberOfMothers") value = object->genNumberOfMothers;
1776 else if(variable == "genMotherId") value = object->genMotherId;
1777 else if(variable == "genMotherCharge") value = object->genMotherCharge;
1778 else if(variable == "genMother0Id") value = object->genMother0Id;
1779 else if(variable == "genMother1Id") value = object->genMother1Id;
1780 else if(variable == "genGrandMother00Id") value = object->genGrandMother00Id;
1781 else if(variable == "genGrandMother01Id") value = object->genGrandMother01Id;
1782 else if(variable == "genGrandMother10Id") value = object->genGrandMother10Id;
1783 else if(variable == "genGrandMother11Id") value = object->genGrandMother11Id;
1784 else if(variable == "isPFMuon") value = object->isPFMuon;
1785 else if(variable == "isGoodMuon_1StationTight") value = object->isGoodMuon_1StationTight;
1786 else if(variable == "isGlobalMuon") value = object->isGlobalMuon;
1787 else if(variable == "isTrackerMuon") value = object->isTrackerMuon;
1788 else if(variable == "isStandAloneMuon") value = object->isStandAloneMuon;
1789 else if(variable == "isGlobalMuonPromptTight") value = object->isGlobalMuonPromptTight;
1790 else if(variable == "numberOfValidMuonHits") value = object->numberOfValidMuonHits;
1791 else if(variable == "numberOfValidTrackerHits") value = object->numberOfValidTrackerHits;
1792 else if(variable == "numberOfLayersWithMeasurement") value = object->numberOfLayersWithMeasurement;
1793 else if(variable == "pixelLayersWithMeasurement") value = object->pixelLayersWithMeasurement;
1794 else if(variable == "numberOfMatches") value = object->numberOfMatches;
1795 else if(variable == "numberOfValidTrackerHitsInnerTrack") value = object->numberOfValidTrackerHitsInnerTrack;
1796 else if(variable == "numberOfValidPixelHits") value = object->numberOfValidPixelHits;
1797 else if(variable == "numberOfMatchedStations") value = object->numberOfMatchedStations;
1798 else if(variable == "time_ndof") value = object->time_ndof;
1799
1800 //user-defined variables
1801 else if(variable == "looseID") {
1802 value = object->pt > 10 &&
1803 (object->isGlobalMuon > 0 ||
1804 object->isTrackerMuon > 0);
1805 }
1806 else if(variable == "correctedD0VertexErr") value = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1807 else if(variable == "correctedD0VertexSig") value = object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1808 else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1809 else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1810 else if(variable == "relPFdBetaIsoPseudo") value = (object->pfIsoR04SumChargedHadronPt + object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt) / object->pt;
1811 else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1812 else if(variable == "metMT") {
1813 if (const BNmet *met = chosenMET ())
1814 {
1815 double dPhi = deltaPhi (object->phi, met->phi);
1816 value = sqrt (2 * object->pt * met->pt * (1 - cos (dPhi)));
1817 }
1818 else
1819 value = -999;
1820 }
1821
1822
1823
1824 else if(variable == "correctedD0VertexInEBPlus"){
1825 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1826 else value = -999;
1827 }
1828 else if(variable == "correctedD0VertexOutEBPlus"){
1829 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1830 else value = -999;
1831 }
1832 else if(variable == "correctedD0VertexEEPlus"){
1833 if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex;
1834 else value = -999;
1835 }
1836
1837 else if(variable == "correctedD0BeamspotInEBPlus"){
1838 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1839 else value = -999;
1840 }
1841 else if(variable == "correctedD0BeamspotOutEBPlus"){
1842 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1843 else value = -999;
1844 }
1845 else if(variable == "correctedD0BeamspotEEPlus"){
1846 if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0;
1847 else value = -999;
1848 }
1849
1850 else if(variable == "correctedD0VertexInEBMinus"){
1851 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1852 else value = -999;
1853 }
1854 else if(variable == "correctedD0VertexOutEBMinus"){
1855 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1856 else value = -999;
1857 }
1858 else if(variable == "correctedD0VertexEEMinus"){
1859 if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex;
1860 else value = -999;
1861 }
1862
1863 else if(variable == "correctedD0BeamspotInEBMinus"){
1864 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1865 else value = -999;
1866 }
1867 else if(variable == "correctedD0BeamspotOutEBMinus"){
1868 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1869 else value = -999;
1870 }
1871 else if(variable == "correctedD0BeamspotEEMinus"){
1872 if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0;
1873 else value = -999;
1874 }
1875
1876
1877 else if(variable == "correctedD0VertexInEBPositiveCharge"){
1878 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1879 else value = -999;
1880 }
1881 else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1882 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1883 else value = -999;
1884 }
1885 else if(variable == "correctedD0VertexEEPositiveCharge"){
1886 if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex;
1887 else value = -999;
1888 }
1889
1890 else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1891 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1892 else value = -999;
1893 }
1894 else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1895 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1896 else value = -999;
1897 }
1898 else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1899 if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0;
1900 else value = -999;
1901 }
1902
1903 else if(variable == "correctedD0VertexInEBNegativeCharge"){
1904 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1905 else value = -999;
1906 }
1907 else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1908 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1909 else value = -999;
1910 }
1911 else if(variable == "correctedD0VertexEENegativeCharge"){
1912 if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex;
1913 else value = -999;
1914 }
1915
1916 else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1917 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1918 else value = -999;
1919 }
1920 else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1921 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1922 else value = -999;
1923 }
1924 else if(variable == "correctedD0BeamspotEENegativeCharge"){
1925 if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0;
1926 else value = -999;
1927 }
1928
1929
1930 else if(variable == "tightID") {
1931 value = object->isGlobalMuon > 0 \
1932 && object->isPFMuon > 0 \
1933 && object->normalizedChi2 < 10 \
1934 && object->numberOfValidMuonHits > 0 \
1935 && object->numberOfMatchedStations > 1 \
1936 && fabs(object->correctedD0Vertex) < 0.2 \
1937 && fabs(object->correctedDZ) < 0.5 \
1938 && object->numberOfValidPixelHits > 0 \
1939 && object->numberOfLayersWithMeasurement > 5;
1940 }
1941 else if(variable == "tightIDdisplaced"){
1942 value = object->isGlobalMuon > 0 \
1943 && object->isPFMuon > 0 \
1944 && object->normalizedChi2 < 10 \
1945 && object->numberOfValidMuonHits > 0 \
1946 && object->numberOfMatchedStations > 1 \
1947 && object->numberOfValidPixelHits > 0 \
1948 && object->numberOfLayersWithMeasurement > 5;
1949 }
1950
1951 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1952
1953 else if(variable == "genMatchedPdgId"){
1954 int index = getGenMatchedParticleIndex(object);
1955 if(index == -1) value = 0;
1956 else value = mcparticles->at(index).id;
1957 }
1958
1959 else if(variable == "genMatchedId"){
1960 int index = getGenMatchedParticleIndex(object);
1961 if(index == -1) value = 0;
1962 else value = getPdgIdBinValue(mcparticles->at(index).id);
1963 }
1964 else if(variable == "genMatchedMotherId"){
1965 int index = getGenMatchedParticleIndex(object);
1966 if(index == -1) value = 0;
1967 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1968 }
1969 else if(variable == "genMatchedMotherIdReverse"){
1970 int index = getGenMatchedParticleIndex(object);
1971 if(index == -1) value = 24;
1972 else value = 24 - getPdgIdBinValue(mcparticles->at(index).motherId);
1973 }
1974 else if(variable == "genMatchedGrandmotherId"){
1975 int index = getGenMatchedParticleIndex(object);
1976 if(index == -1) value = 0;
1977 else if(fabs(mcparticles->at(index).motherId) == 15){
1978 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1979 if(motherIndex == -1) value = 0;
1980 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1981 }
1982 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1983 }
1984 else if(variable == "genMatchedGrandmotherIdReverse"){
1985 int index = getGenMatchedParticleIndex(object);
1986 if(index == -1) value = 24;
1987 else if(fabs(mcparticles->at(index).motherId) == 15){
1988 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1989 if(motherIndex == -1) value = 24;
1990 else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1991 }
1992 else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1993 }
1994 else if(variable == "pfMuonsFromVertex"){
1995 double d0Error, dzError;
1996
1997 d0Error = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1998 dzError = hypot (object->tkDZerr, chosenVertex ()->zError);
1999 value = fabs (object->correctedD0Vertex) > 0.2 || fabs (object->correctedDZ) > 0.5
2000 || fabs (object->correctedD0Vertex / d0Error) > 99.0
2001 || fabs (object->correctedDZ / dzError) > 99.0;
2002 value = (object->isStandAloneMuon && !object->isTrackerMuon && !object->isGlobalMuon) || !value;
2003 }
2004
2005
2006
2007 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2008
2009 value = applyFunction(function, value);
2010
2011 return value;
2012 } // end muon valueLookup
2013
2014
2015 //!electron valueLookup
2016 double
2017 OSUAnalysis::valueLookup (const BNelectron* object, string variable, string function, string &stringValue){
2018
2019 double value = 0.0;
2020 if(variable == "energy") value = object->energy;
2021 else if(variable == "et") value = object->et;
2022 else if(variable == "gsfEt") value = object->gsfEt;
2023 else if(variable == "pt") value = object->pt;
2024 else if(variable == "px") value = object->px;
2025 else if(variable == "py") value = object->py;
2026 else if(variable == "pz") value = object->pz;
2027 else if(variable == "phi") value = object->phi;
2028 else if(variable == "eta") value = object->eta;
2029 else if(variable == "theta") value = object->theta;
2030 else if(variable == "pIn") value = object->pIn;
2031 else if(variable == "pOut") value = object->pOut;
2032 else if(variable == "EscOverPin") value = object->EscOverPin;
2033 else if(variable == "EseedOverPout") value = object->EseedOverPout;
2034 else if(variable == "hadOverEm") value = object->hadOverEm;
2035 else if(variable == "trackIso") value = object->trackIso;
2036 else if(variable == "ecalIso") value = object->ecalIso;
2037 else if(variable == "hcalIso") value = object->hcalIso;
2038 else if(variable == "caloIso") value = object->caloIso;
2039 else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
2040 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
2041 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
2042 else if(variable == "hcalIsoDR03depth1") value = object->hcalIsoDR03depth1;
2043 else if(variable == "hcalIsoDR03depth2") value = object->hcalIsoDR03depth2;
2044 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
2045 else if(variable == "trackIsoDR04") value = object->trackIsoDR04;
2046 else if(variable == "ecalIsoDR04") value = object->ecalIsoDR04;
2047 else if(variable == "hcalIsoDR04") value = object->hcalIsoDR04;
2048 else if(variable == "hcalIsoDR04depth1") value = object->hcalIsoDR04depth1;
2049 else if(variable == "hcalIsoDR04depth2") value = object->hcalIsoDR04depth2;
2050 else if(variable == "caloIsoDR04") value = object->caloIsoDR04;
2051 else if(variable == "fbrem") value = object->fbrem;
2052 else if(variable == "absInvEMinusInvPin") value = object->absInvEMinusInvPin;
2053 else if(variable == "delPhiIn") value = object->delPhiIn;
2054 else if(variable == "delEtaIn") value = object->delEtaIn;
2055 else if(variable == "genET") value = object->genET;
2056 else if(variable == "genPT") value = object->genPT;
2057 else if(variable == "genPhi") value = object->genPhi;
2058 else if(variable == "genEta") value = object->genEta;
2059 else if(variable == "genMotherET") value = object->genMotherET;
2060 else if(variable == "genMotherPT") value = object->genMotherPT;
2061 else if(variable == "genMotherPhi") value = object->genMotherPhi;
2062 else if(variable == "genMotherEta") value = object->genMotherEta;
2063 else if(variable == "vx") value = object->vx;
2064 else if(variable == "vy") value = object->vy;
2065 else if(variable == "vz") value = object->vz;
2066 else if(variable == "scEnergy") value = object->scEnergy;
2067 else if(variable == "scRawEnergy") value = object->scRawEnergy;
2068 else if(variable == "scSigmaEtaEta") value = object->scSigmaEtaEta;
2069 else if(variable == "scSigmaIEtaIEta") value = object->scSigmaIEtaIEta;
2070 else if(variable == "scE1x5") value = object->scE1x5;
2071 else if(variable == "scE2x5Max") value = object->scE2x5Max;
2072 else if(variable == "scE5x5") value = object->scE5x5;
2073 else if(variable == "scEt") value = object->scEt;
2074 else if(variable == "scEta") value = object->scEta;
2075 else if(variable == "scPhi") value = object->scPhi;
2076 else if(variable == "scZ") value = object->scZ;
2077 else if(variable == "tkNormChi2") value = object->tkNormChi2;
2078 else if(variable == "tkPT") value = object->tkPT;
2079 else if(variable == "tkEta") value = object->tkEta;
2080 else if(variable == "tkPhi") value = object->tkPhi;
2081 else if(variable == "tkDZ") value = object->tkDZ;
2082 else if(variable == "tkD0") value = object->tkD0;
2083 else if(variable == "tkD0bs") value = object->tkD0bs;
2084 else if(variable == "tkD0err") value = object->tkD0err;
2085 else if(variable == "mva") value = object->mva;
2086 else if(variable == "mvaTrigV0") value = object->mvaTrigV0;
2087 else if(variable == "mvaNonTrigV0") value = object->mvaNonTrigV0;
2088 else if(variable == "dist") value = object->dist;
2089 else if(variable == "dcot") value = object->dcot;
2090 else if(variable == "convradius") value = object->convradius;
2091 else if(variable == "convPointX") value = object->convPointX;
2092 else if(variable == "convPointY") value = object->convPointY;
2093 else if(variable == "convPointZ") value = object->convPointZ;
2094 else if(variable == "eMax") value = object->eMax;
2095 else if(variable == "eLeft") value = object->eLeft;
2096 else if(variable == "eRight") value = object->eRight;
2097 else if(variable == "eTop") value = object->eTop;
2098 else if(variable == "eBottom") value = object->eBottom;
2099 else if(variable == "e3x3") value = object->e3x3;
2100 else if(variable == "swissCross") value = object->swissCross;
2101 else if(variable == "seedEnergy") value = object->seedEnergy;
2102 else if(variable == "seedTime") value = object->seedTime;
2103 else if(variable == "swissCrossNoI85") value = object->swissCrossNoI85;
2104 else if(variable == "swissCrossI85") value = object->swissCrossI85;
2105 else if(variable == "E2overE9NoI85") value = object->E2overE9NoI85;
2106 else if(variable == "E2overE9I85") value = object->E2overE9I85;
2107 else if(variable == "correctedD0") value = object->correctedD0;
2108 else if(variable == "correctedD0Vertex") value = object->correctedD0Vertex;
2109 else if(variable == "correctedDZ") value = object->correctedDZ;
2110 else if(variable == "particleIso") value = object->particleIso;
2111 else if(variable == "chargedHadronIso") value = object->chargedHadronIso;
2112 else if(variable == "neutralHadronIso") value = object->neutralHadronIso;
2113 else if(variable == "photonIso") value = object->photonIso;
2114 else if(variable == "puChargedHadronIso") value = object->puChargedHadronIso;
2115 else if(variable == "chargedHadronIsoDR03") value = object->chargedHadronIsoDR03;
2116 else if(variable == "neutralHadronIsoDR03") value = object->neutralHadronIsoDR03;
2117 else if(variable == "photonIsoDR03") value = object->photonIsoDR03;
2118 else if(variable == "puChargedHadronIsoDR03") value = object->puChargedHadronIsoDR03;
2119 else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
2120 else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
2121 else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
2122 else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
2123 else if(variable == "rhoPrime") value = object->rhoPrime;
2124 else if(variable == "AEffDr03") value = object->AEffDr03;
2125 else if(variable == "AEffDr04") value = object->AEffDr04;
2126 else if(variable == "IP") value = object->IP;
2127 else if(variable == "IPError") value = object->IPError;
2128 else if(variable == "charge") value = object->charge;
2129 else if(variable == "classification") value = object->classification;
2130 else if(variable == "genId") value = object->genId;
2131 else if(variable == "genCharge") value = object->genCharge;
2132 else if(variable == "genNumberOfMothers") value = object->genNumberOfMothers;
2133 else if(variable == "genMotherId") value = object->genMotherId;
2134 else if(variable == "genMotherCharge") value = object->genMotherCharge;
2135 else if(variable == "genMother0Id") value = object->genMother0Id;
2136 else if(variable == "genMother1Id") value = object->genMother1Id;
2137 else if(variable == "genGrandMother00Id") value = object->genGrandMother00Id;
2138 else if(variable == "genGrandMother01Id") value = object->genGrandMother01Id;
2139 else if(variable == "genGrandMother10Id") value = object->genGrandMother10Id;
2140 else if(variable == "genGrandMother11Id") value = object->genGrandMother11Id;
2141 else if(variable == "numClusters") value = object->numClusters;
2142 else if(variable == "tkNumValidHits") value = object->tkNumValidHits;
2143 else if(variable == "tkCharge") value = object->tkCharge;
2144 else if(variable == "gsfCharge") value = object->gsfCharge;
2145 else if(variable == "isEB") value = object->isEB;
2146 else if(variable == "isEE") value = object->isEE;
2147 else if(variable == "isGap") value = object->isGap;
2148 else if(variable == "isEBEEGap") value = object->isEBEEGap;
2149 else if(variable == "isEBGap") value = object->isEBGap;
2150 else if(variable == "isEEGap") value = object->isEEGap;
2151 else if(variable == "isEcalDriven") value = object->isEcalDriven;
2152 else if(variable == "isTrackerDriven") value = object->isTrackerDriven;
2153 else if(variable == "numberOfLostHits") value = object->numberOfLostHits;
2154 else if(variable == "numberOfExpectedInnerHits") value = object->numberOfExpectedInnerHits;
2155 else if(variable == "numberOfValidPixelHits") value = object->numberOfValidPixelHits;
2156 else if(variable == "numberOfValidPixelBarrelHits") value = object->numberOfValidPixelBarrelHits;
2157 else if(variable == "numberOfValidPixelEndcapHits") value = object->numberOfValidPixelEndcapHits;
2158 else if(variable == "isHEEP") value = object->isHEEP;
2159 else if(variable == "isHEEPnoEt") value = object->isHEEPnoEt;
2160 else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
2161 else if(variable == "eidRobustHighEnergy") value = object->eidRobustHighEnergy;
2162 else if(variable == "eidRobustLoose") value = object->eidRobustLoose;
2163 else if(variable == "eidRobustTight") value = object->eidRobustTight;
2164 else if(variable == "eidLoose") value = object->eidLoose;
2165 else if(variable == "eidTight") value = object->eidTight;
2166 else if(variable == "eidVeryLooseMC") value = object->eidVeryLooseMC;
2167 else if(variable == "eidLooseMC") value = object->eidLooseMC;
2168 else if(variable == "eidMediumMC") value = object->eidMediumMC;
2169 else if(variable == "eidTightMC") value = object->eidTightMC;
2170 else if(variable == "eidSuperTightMC") value = object->eidSuperTightMC;
2171 else if(variable == "eidHyperTight1MC") value = object->eidHyperTight1MC;
2172 else if(variable == "eidHyperTight2MC") value = object->eidHyperTight2MC;
2173 else if(variable == "eidHyperTight3MC") value = object->eidHyperTight3MC;
2174 else if(variable == "eidHyperTight4MC") value = object->eidHyperTight4MC;
2175 else if(variable == "passConvVeto") value = object->passConvVeto;
2176
2177 //user-defined variables
2178 else if(variable == "correctedD0VertexErr") value = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
2179 else if(variable == "correctedD0VertexSig") value = object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
2180 else if(variable == "detIso") value = (object->trackIso) / object->pt;
2181 else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
2182 else if(variable == "relPFrhoIsoEB") value = object->isEB ? ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt : -999;
2183 else if(variable == "relPFrhoIsoEE") value = object->isEE ? ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt : -999;
2184 else if(variable == "metMT") {
2185 if (const BNmet *met = chosenMET ())
2186 {
2187 double dPhi = deltaPhi (object->phi, met->phi);
2188 value = sqrt (2 * object->pt * met->pt * (1 - cos (dPhi)));
2189 }
2190 else
2191 value = -999;
2192 }
2193
2194 else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
2195 if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
2196 else value = -999;
2197 }
2198 else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){
2199 if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex;
2200 else value = -999;
2201 }
2202
2203 else if(variable == "correctedD0VertexInEBPlus"){
2204 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
2205 else value = -999;
2206 }
2207 else if(variable == "correctedD0VertexOutEBPlus"){
2208 if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
2209 else value = -999;
2210 }
2211 else if(variable == "correctedD0VertexEEPlus"){
2212 if(object->isEE && object->eta > 0) value = object->correctedD0Vertex;
2213 else value = -999;
2214 }
2215
2216 else if(variable == "correctedD0BeamspotInEBPlus"){
2217 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
2218 else value = -999;
2219 }
2220 else if(variable == "correctedD0BeamspotOutEBPlus"){
2221 if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
2222 else value = -999;
2223 }
2224 else if(variable == "correctedD0BeamspotEEPlus"){
2225 if(object->isEE && object->eta > 0) value = object->correctedD0;
2226 else value = -999;
2227 }
2228
2229 else if(variable == "correctedD0VertexInEBMinus"){
2230 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
2231 else value = -999;
2232 }
2233 else if(variable == "correctedD0VertexOutEBMinus"){
2234 if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
2235 else value = -999;
2236 }
2237 else if(variable == "correctedD0VertexEEMinus"){
2238 if(object->isEE && object->eta < 0) value = object->correctedD0Vertex;
2239 else value = -999;
2240 }
2241
2242 else if(variable == "correctedD0BeamspotInEBMinus"){
2243 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
2244 else value = -999;
2245 }
2246 else if(variable == "correctedD0BeamspotOutEBMinus"){
2247 if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
2248 else value = -999;
2249 }
2250 else if(variable == "correctedD0BeamspotEEMinus"){
2251 if(object->isEE && object->eta < 0) value = object->correctedD0;
2252 else value = -999;
2253 }
2254
2255 else if(variable == "looseID"){
2256 if (object->isEB)
2257 {
2258 value = fabs(object->delEtaIn) < 0.007 \
2259 && fabs (object->delPhiIn) < 0.15 \
2260 && object->scSigmaIEtaIEta < 0.01 \
2261 && object->hadOverEm < 0.12 \
2262 && fabs (object->correctedD0Vertex) < 0.02 \
2263 && fabs (object->correctedDZ) < 0.2 \
2264 && object->absInvEMinusInvPin < 0.05 \
2265 && object->passConvVeto;
2266 }
2267 else
2268 {
2269 value = fabs(object->delEtaIn) < 0.009 \
2270 && fabs (object->delPhiIn) < 0.10 \
2271 && object->scSigmaIEtaIEta < 0.03 \
2272 && object->hadOverEm < 0.10 \
2273 && fabs (object->correctedD0Vertex) < 0.02 \
2274 && fabs (object->correctedDZ) < 0.2 \
2275 && object->absInvEMinusInvPin < 0.05 \
2276 && object->passConvVeto;
2277 }
2278 }
2279
2280 else if(variable == "tightID"){
2281 if (object->isEB)
2282 {
2283 value = fabs(object->delEtaIn) < 0.004 \
2284 && fabs (object->delPhiIn) < 0.03 \
2285 && object->scSigmaIEtaIEta < 0.01 \
2286 && object->hadOverEm < 0.12 \
2287 && fabs (object->correctedD0Vertex) < 0.02 \
2288 && fabs (object->correctedDZ) < 0.1 \
2289 && object->absInvEMinusInvPin < 0.05 \
2290 && object->passConvVeto;
2291 }
2292 else
2293 {
2294 value = fabs(object->delEtaIn) < 0.005 \
2295 && fabs (object->delPhiIn) < 0.02 \
2296 && object->scSigmaIEtaIEta < 0.03 \
2297 && object->hadOverEm < 0.10 \
2298 && fabs (object->correctedD0Vertex) < 0.02 \
2299 && fabs (object->correctedDZ) < 0.1 \
2300 && object->absInvEMinusInvPin < 0.05 \
2301 && object->passConvVeto;
2302 }
2303 }
2304
2305 else if(variable == "looseID_MVA"){
2306 value = object->pt > 10
2307 && object->mvaNonTrigV0 > 0;
2308 }
2309 else if(variable == "correctedD0VertexInEBPositiveCharge"){
2310 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
2311 else value = -999;
2312 }
2313 else if(variable == "correctedD0VertexOutEBPositiveCharge"){
2314 if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
2315 else value = -999;
2316 }
2317 else if(variable == "correctedD0VertexEEPositiveCharge"){
2318 if(object->isEE && object->charge > 0) value = object->correctedD0Vertex;
2319 else value = -999;
2320 }
2321
2322 else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
2323 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
2324 else value = -999;
2325 }
2326 else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
2327 if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
2328 else value = -999;
2329 }
2330 else if(variable == "correctedD0BeamspotEEPositiveCharge"){
2331 if(object->isEE && object->charge > 0) value = object->correctedD0;
2332 else value = -999;
2333 }
2334
2335 else if(variable == "correctedD0VertexInEBNegativeCharge"){
2336 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
2337 else value = -999;
2338 }
2339 else if(variable == "correctedD0VertexOutEBNegativeCharge"){
2340 if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
2341 else value = -999;
2342 }
2343 else if(variable == "correctedD0VertexEENegativeCharge"){
2344 if(object->isEE && object->charge < 0) value = object->correctedD0Vertex;
2345 else value = -999;
2346 }
2347
2348 else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
2349 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
2350 else value = -999;
2351 }
2352 else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
2353 if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
2354 else value = -999;
2355 }
2356 else if(variable == "correctedD0BeamspotEENegativeCharge"){
2357 if(object->isEE && object->charge < 0) value = object->correctedD0;
2358 else value = -999;
2359 }
2360
2361
2362 else if(variable == "tightIDdisplaced"){
2363 if (object->isEB)
2364 {
2365 value = fabs(object->delEtaIn) < 0.004 \
2366 && fabs (object->delPhiIn) < 0.03 \
2367 && object->scSigmaIEtaIEta < 0.01 \
2368 && object->hadOverEm < 0.12 \
2369 && object->absInvEMinusInvPin < 0.05;
2370 }
2371 else
2372 {
2373 value = fabs (object->delEtaIn) < 0.005 \
2374 && fabs (object->delPhiIn) < 0.02 \
2375 && object->scSigmaIEtaIEta < 0.03 \
2376 && object->hadOverEm < 0.10 \
2377 && object->absInvEMinusInvPin < 0.05;
2378 }
2379 }
2380
2381
2382 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2383
2384 else if(variable == "genMatchedPdgId"){
2385 int index = getGenMatchedParticleIndex(object);
2386 if(index == -1) value = 0;
2387 else value = mcparticles->at(index).id;
2388 }
2389
2390
2391 else if(variable == "genMatchedId"){
2392 int index = getGenMatchedParticleIndex(object);
2393 if(index == -1) value = 0;
2394 else value = getPdgIdBinValue(mcparticles->at(index).id);
2395 }
2396 else if(variable == "genMatchedMotherId"){
2397 int index = getGenMatchedParticleIndex(object);
2398 if(index == -1) value = 0;
2399 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2400 }
2401 else if(variable == "genMatchedMotherIdReverse"){
2402 int index = getGenMatchedParticleIndex(object);
2403 if(index == -1) value = 24;
2404 else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2405 }
2406 else if(variable == "genMatchedGrandmotherId"){
2407 int index = getGenMatchedParticleIndex(object);
2408 if(index == -1) value = 0;
2409 else if(fabs(mcparticles->at(index).motherId) == 15){
2410 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2411 if(motherIndex == -1) value = 0;
2412 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2413 }
2414 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2415 }
2416 else if(variable == "genMatchedGrandmotherIdReverse"){
2417 int index = getGenMatchedParticleIndex(object);
2418 if(index == -1) value = 24;
2419 else if(fabs(mcparticles->at(index).motherId) == 15){
2420 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2421 if(motherIndex == -1) value = 24;
2422 else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2423 }
2424 else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2425 }
2426 else if(variable == "pfElectronsFromVertex"){
2427 double d0Error, dzError;
2428
2429 d0Error = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
2430 dzError = hypot (object->tkDZerr, chosenVertex ()->zError);
2431 value = fabs (object->correctedD0Vertex) > 0.2 || fabs (object->correctedDZ) > 0.5
2432 || fabs (object->correctedD0Vertex / d0Error) > 99.0
2433 || fabs (object->correctedDZ / dzError) > 99.0;
2434 value = !value;
2435 }
2436
2437
2438
2439 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2440
2441 value = applyFunction(function, value);
2442
2443 return value;
2444 } // end electron valueLookup
2445
2446
2447 //!event valueLookup
2448 double
2449 OSUAnalysis::valueLookup (const BNevent* object, string variable, string function, string &stringValue){
2450
2451 double value = 0.0;
2452
2453 if(variable == "weight") value = object->weight;
2454 else if(variable == "pthat") value = object->pthat;
2455 else if(variable == "qScale") value = object->qScale;
2456 else if(variable == "alphaQCD") value = object->alphaQCD;
2457 else if(variable == "Ht") value = getHt(jets.product());
2458 else if(variable == "alphaQED") value = object->alphaQED;
2459 else if(variable == "scalePDF") value = object->scalePDF;
2460 else if(variable == "x1") value = object->x1;
2461 else if(variable == "x2") value = object->x2;
2462 else if(variable == "xPDF1") value = object->xPDF1;
2463 else if(variable == "xPDF2") value = object->xPDF2;
2464 else if(variable == "BSx") value = object->BSx;
2465 else if(variable == "BSy") value = object->BSy;
2466 else if(variable == "BSz") value = object->BSz;
2467 else if(variable == "PVx") value = object->PVx;
2468 else if(variable == "PVy") value = object->PVy;
2469 else if(variable == "PVz") value = object->PVz;
2470 else if(variable == "bField") value = object->bField;
2471 else if(variable == "instLumi") value = object->instLumi;
2472 else if(variable == "bxLumi") value = object->bxLumi;
2473 else if(variable == "FilterOutScrapingFraction") value = object->FilterOutScrapingFraction;
2474 else if(variable == "sumNVtx") value = object->sumNVtx;
2475 else if(variable == "sumTrueNVtx") value = object->sumTrueNVtx;
2476 else if(variable == "nm1_true") value = object->nm1_true;
2477 else if(variable == "n0_true") value = object->n0_true;
2478 else if(variable == "np1_true") value = object->np1_true;
2479 else if(variable == "numTruePV") value = object->numTruePV;
2480 else if(variable == "Q2ScaleUpWgt") value = object->Q2ScaleUpWgt;
2481 else if(variable == "Q2ScaleDownWgt") value = object->Q2ScaleDownWgt;
2482 else if(variable == "rho_kt6PFJets") value = object->rho_kt6PFJets;
2483 else if(variable == "rho_kt6PFJetsCentralChargedPileUp") value = object->rho_kt6PFJetsCentralChargedPileUp;
2484 else if(variable == "rho_kt6PFJetsCentralNeutral") value = object->rho_kt6PFJetsCentralNeutral;
2485 else if(variable == "rho_kt6PFJetsCentralNeutralTight") value = object->rho_kt6PFJetsCentralNeutralTight;
2486 else if(variable == "run") value = object->run;
2487 else if(variable == "lumi") value = object->lumi;
2488 else if(variable == "sample") value = object->sample;
2489 else if(variable == "numPV") value = object->numPV;
2490 else if(variable == "W0decay") value = object->W0decay;
2491 else if(variable == "W1decay") value = object->W1decay;
2492 else if(variable == "Z0decay") value = object->Z0decay;
2493 else if(variable == "Z1decay") value = object->Z1decay;
2494 else if(variable == "H0decay") value = object->H0decay;
2495 else if(variable == "H1decay") value = object->H1decay;
2496 else if(variable == "hcalnoiseLoose") value = object->hcalnoiseLoose;
2497 else if(variable == "hcalnoiseTight") value = object->hcalnoiseTight;
2498 else if(variable == "GoodVertex") value = object->GoodVertex;
2499 else if(variable == "FilterOutScraping") value = object->FilterOutScraping;
2500 else if(variable == "HBHENoiseFilter") value = object->HBHENoiseFilter;
2501 else if(variable == "CSCLooseHaloId") value = object->CSCLooseHaloId;
2502 else if(variable == "CSCTightHaloId") value = object->CSCTightHaloId;
2503 else if(variable == "EcalLooseHaloId") value = object->EcalLooseHaloId;
2504 else if(variable == "EcalTightHaloId") value = object->EcalTightHaloId;
2505 else if(variable == "HcalLooseHaloId") value = object->HcalLooseHaloId;
2506 else if(variable == "HcalTightHaloId") value = object->HcalTightHaloId;
2507 else if(variable == "GlobalLooseHaloId") value = object->GlobalLooseHaloId;
2508 else if(variable == "GlobalTightHaloId") value = object->GlobalTightHaloId;
2509 else if(variable == "LooseId") value = object->LooseId;
2510 else if(variable == "TightId") value = object->TightId;
2511 else if(variable == "numGenPV") value = object->numGenPV;
2512 else if(variable == "nm1") value = object->nm1;
2513 else if(variable == "n0") value = object->n0;
2514 else if(variable == "np1") value = object->np1;
2515 else if(variable == "id1") value = object->id1;
2516 else if(variable == "id2") value = object->id2;
2517 else if(variable == "evt") value = object->evt;
2518 else if(variable == "puScaleFactor"){
2519 if(doPileupReweighting_ && datasetType_ != "data")
2520 value = puWeight_->at (events->at (0).numTruePV);
2521 else
2522 value = 1.0;
2523 }
2524 else if(variable == "muonScaleFactor") value = muonScaleFactor_;
2525 else if(variable == "electronScaleFactor") value = electronScaleFactor_;
2526 else if(variable == "stopCTauScaleFactor") value = stopCTauScaleFactor_;
2527 else if(variable == "bTagScaleFactor") value = bTagScaleFactor_;
2528 else if(variable == "ht") value = chosenHT ();
2529 else if(variable == "leadMuPairInvMass"){
2530 pair<const BNmuon *, const BNmuon *> muPair = leadMuonPair ();
2531 TLorentzVector p0 (muPair.first->px, muPair.first->py, muPair.first->pz, muPair.first->energy),
2532 p1 (muPair.second->px, muPair.second->py, muPair.second->pz, muPair.second->energy);
2533 value = (p0 + p1).M ();
2534 }
2535 else if(variable == "leadMuPairPt"){
2536 pair<const BNmuon *, const BNmuon *> muPair = leadMuonPair ();
2537 TVector2 pt0 (muPair.first->px, muPair.first->py),
2538 pt1 (muPair.second->px, muPair.second->py);
2539 pt0 += pt1;
2540 value = pt0.Mod ();
2541 }
2542 else if(variable == "leadElPairInvMass"){
2543 pair<const BNelectron *, const BNelectron *> muPair = leadElectronPair ();
2544 TLorentzVector p0 (muPair.first->px, muPair.first->py, muPair.first->pz, muPair.first->energy),
2545 p1 (muPair.second->px, muPair.second->py, muPair.second->pz, muPair.second->energy);
2546 value = (p0 + p1).M ();
2547 }
2548 else if(variable == "leadElPairPt"){
2549 pair<const BNelectron *, const BNelectron *> muPair = leadElectronPair ();
2550 TVector2 pt0 (muPair.first->px, muPair.first->py),
2551 pt1 (muPair.second->px, muPair.second->py);
2552 pt0 += pt1;
2553 value = pt0.Mod ();
2554 }
2555 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2556
2557 value = applyFunction(function, value);
2558
2559 return value;
2560 } // end event valueLookup
2561
2562
2563 //!tau valueLookup
2564 double
2565 OSUAnalysis::valueLookup (const BNtau* object, string variable, string function, string &stringValue){
2566
2567 double value = 0.0;
2568
2569 if(variable == "px") value = object->px;
2570 else if(variable == "py") value = object->py;
2571 else if(variable == "pz") value = object->pz;
2572 else if(variable == "energy") value = object->energy;
2573 else if(variable == "et") value = object->et;
2574 else if(variable == "pt") value = object->pt;
2575 else if(variable == "eta") value = object->eta;
2576 else if(variable == "phi") value = object->phi;
2577 else if(variable == "emFraction") value = object->emFraction;
2578 else if(variable == "leadingTrackPt") value = object->leadingTrackPt;
2579 else if(variable == "leadingTrackIpVtdxy") value = object->leadingTrackIpVtdxy;
2580 else if(variable == "leadingTrackIpVtdz") value = object->leadingTrackIpVtdz;
2581 else if(variable == "leadingTrackIpVtdxyError") value = object->leadingTrackIpVtdxyError;
2582 else if(variable == "leadingTrackIpVtdzError") value = object->leadingTrackIpVtdzError;
2583 else if(variable == "leadingTrackVx") value = object->leadingTrackVx;
2584 else if(variable == "leadingTrackVy") value = object->leadingTrackVy;
2585 else if(variable == "leadingTrackVz") value = object->leadingTrackVz;
2586 else if(variable == "leadingTrackValidHits") value = object->leadingTrackValidHits;
2587 else if(variable == "leadingTrackNormChiSqrd") value = object->leadingTrackNormChiSqrd;
2588 else if(variable == "numProngs") value = object->numProngs;
2589 else if(variable == "numSignalGammas") value = object->numSignalGammas;
2590 else if(variable == "numSignalNeutrals") value = object->numSignalNeutrals;
2591 else if(variable == "numSignalPiZeros") value = object->numSignalPiZeros;
2592 else if(variable == "decayMode") value = object->decayMode;
2593 else if(variable == "charge") value = object->charge;
2594 else if(variable == "inTheCracks") value = object->inTheCracks;
2595 else if(variable == "HPSagainstElectronLoose") value = object->HPSagainstElectronLoose;
2596 else if(variable == "HPSagainstElectronMVA") value = object->HPSagainstElectronMVA;
2597 else if(variable == "HPSagainstElectronMedium") value = object->HPSagainstElectronMedium;
2598 else if(variable == "HPSagainstElectronTight") value = object->HPSagainstElectronTight;
2599 else if(variable == "HPSagainstMuonLoose") value = object->HPSagainstMuonLoose;
2600 else if(variable == "HPSagainstMuonMedium") value = object->HPSagainstMuonMedium;
2601 else if(variable == "HPSagainstMuonTight") value = object->HPSagainstMuonTight;
2602 else if(variable == "HPSbyLooseCombinedIsolationDeltaBetaCorr") value = object->HPSbyLooseCombinedIsolationDeltaBetaCorr;
2603 else if(variable == "HPSbyMediumCombinedIsolationDeltaBetaCorr") value = object->HPSbyMediumCombinedIsolationDeltaBetaCorr;
2604 else if(variable == "HPSbyTightCombinedIsolationDeltaBetaCorr") value = object->HPSbyTightCombinedIsolationDeltaBetaCorr;
2605 else if(variable == "HPSbyVLooseCombinedIsolationDeltaBetaCorr") value = object->HPSbyVLooseCombinedIsolationDeltaBetaCorr;
2606 else if(variable == "HPSdecayModeFinding") value = object->HPSdecayModeFinding;
2607 else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
2608
2609 else if (variable == "looseHadronicID") {
2610 value = object->pt > 10
2611 && object->eta < 2.3
2612 && object->HPSbyLooseCombinedIsolationDeltaBetaCorr > 0
2613 && object->HPSdecayModeFinding > 0
2614 && object->HPSagainstElectronLoose > 0
2615 && object->HPSagainstMuonTight > 0;
2616 }
2617
2618 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2619
2620 else if(variable == "genMatchedPdgId"){
2621 int index = getGenMatchedParticleIndex(object);
2622 if(index == -1) value = 0;
2623 else value = mcparticles->at(index).id;
2624 }
2625
2626 else if(variable == "genMatchedId"){
2627 int index = getGenMatchedParticleIndex(object);
2628 if(index == -1) value = 0;
2629 else value = getPdgIdBinValue(mcparticles->at(index).id);
2630 }
2631 else if(variable == "genMatchedMotherId"){
2632 int index = getGenMatchedParticleIndex(object);
2633 if(index == -1) value = 0;
2634 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2635 }
2636 else if(variable == "genMatchedMotherIdReverse"){
2637 int index = getGenMatchedParticleIndex(object);
2638 if(index == -1) value = 24;
2639 else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2640 }
2641 else if(variable == "genMatchedGrandmotherId"){
2642 int index = getGenMatchedParticleIndex(object);
2643 if(index == -1) value = 0;
2644 else if(fabs(mcparticles->at(index).motherId) == 15){
2645 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2646 if(motherIndex == -1) value = 0;
2647 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2648 }
2649 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2650 }
2651 else if(variable == "genMatchedGrandmotherIdReverse"){
2652 int index = getGenMatchedParticleIndex(object);
2653 if(index == -1) value = 24;
2654 else if(fabs(mcparticles->at(index).motherId) == 15){
2655 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2656 if(motherIndex == -1) value = 24;
2657 else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2658 }
2659 else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2660 }
2661
2662
2663 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2664
2665 value = applyFunction(function, value);
2666
2667 return value;
2668 } // end tau valueLookup
2669
2670
2671 //!met valueLookup
2672 double
2673 OSUAnalysis::valueLookup (const BNmet* object, string variable, string function, string &stringValue){
2674
2675 double value = 0.0;
2676
2677 if(variable == "et") value = object->et;
2678 else if(variable == "pt") value = object->pt;
2679 else if(variable == "px") value = object->px;
2680 else if(variable == "py") value = object->py;
2681 else if(variable == "phi") value = object->phi;
2682 else if(variable == "Upt") value = object->Upt;
2683 else if(variable == "Uphi") value = object->Uphi;
2684 else if(variable == "NeutralEMFraction") value = object->NeutralEMFraction;
2685 else if(variable == "NeutralHadEtFraction") value = object->NeutralHadEtFraction;
2686 else if(variable == "ChargedEMEtFraction") value = object->ChargedEMEtFraction;
2687 else if(variable == "ChargedHadEtFraction") value = object->ChargedHadEtFraction;
2688 else if(variable == "MuonEtFraction") value = object->MuonEtFraction;
2689 else if(variable == "Type6EtFraction") value = object->Type6EtFraction;
2690 else if(variable == "Type7EtFraction") value = object->Type7EtFraction;
2691 else if(variable == "genPT") value = object->genPT;
2692 else if(variable == "genPhi") value = object->genPhi;
2693 else if(variable == "muonCorEx") value = object->muonCorEx;
2694 else if(variable == "muonCorEy") value = object->muonCorEy;
2695 else if(variable == "jet20CorEx") value = object->jet20CorEx;
2696 else if(variable == "jet20CorEy") value = object->jet20CorEy;
2697 else if(variable == "jet1CorEx") value = object->jet1CorEx;
2698 else if(variable == "jet1CorEy") value = object->jet1CorEy;
2699 else if(variable == "sumET") value = object->sumET;
2700 else if(variable == "corSumET") value = object->corSumET;
2701 else if(variable == "mEtSig") value = object->mEtSig;
2702 else if(variable == "metSignificance") value = object->metSignificance;
2703 else if(variable == "significance") value = object->significance;
2704 else if(variable == "sigmaX2") value = object->sigmaX2;
2705 else if(variable == "sigmaY2") value = object->sigmaY2;
2706 else if(variable == "sigmaXY") value = object->sigmaXY;
2707 else if(variable == "sigmaYX") value = object->sigmaYX;
2708 else if(variable == "maxEtInEmTowers") value = object->maxEtInEmTowers;
2709 else if(variable == "emEtFraction") value = object->emEtFraction;
2710 else if(variable == "emEtInEB") value = object->emEtInEB;
2711 else if(variable == "emEtInEE") value = object->emEtInEE;
2712 else if(variable == "emEtInHF") value = object->emEtInHF;
2713 else if(variable == "maxEtInHadTowers") value = object->maxEtInHadTowers;
2714 else if(variable == "hadEtFraction") value = object->hadEtFraction;
2715 else if(variable == "hadEtInHB") value = object->hadEtInHB;
2716 else if(variable == "hadEtInHE") value = object->hadEtInHE;
2717 else if(variable == "hadEtInHF") value = object->hadEtInHF;
2718 else if(variable == "hadEtInHO") value = object->hadEtInHO;
2719 else if(variable == "UDeltaPx") value = object->UDeltaPx;
2720 else if(variable == "UDeltaPy") value = object->UDeltaPy;
2721 else if(variable == "UDeltaP") value = object->UDeltaP;
2722 else if(variable == "Uscale") value = object->Uscale;
2723 else if(variable == "type2corPx") value = object->type2corPx;
2724 else if(variable == "type2corPy") value = object->type2corPy;
2725 else if(variable == "T2pt") value = object->T2pt;
2726 else if(variable == "T2px") value = object->T2px;
2727 else if(variable == "T2py") value = object->T2py;
2728 else if(variable == "T2phi") value = object->T2phi;
2729 else if(variable == "T2sumET") value = object->T2sumET;
2730 else if(variable == "pfT1jet1pt") value = object->pfT1jet1pt;
2731 else if(variable == "pfT1jet1phi") value = object->pfT1jet1phi;
2732 else if(variable == "pfT1jet6pt") value = object->pfT1jet6pt;
2733 else if(variable == "pfT1jet6phi") value = object->pfT1jet6phi;
2734 else if(variable == "pfT1jet10pt") value = object->pfT1jet10pt;
2735 else if(variable == "pfT1jet10phi") value = object->pfT1jet10phi;
2736
2737 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2738
2739 value = applyFunction(function, value);
2740
2741 return value;
2742 } // end met valueLookup
2743
2744
2745 //!track valueLookup
2746 double
2747 OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function, string &stringValue){
2748
2749 double value = 0.0;
2750 double pMag = sqrt(object->pt * object->pt +
2751 object->pz * object->pz);
2752
2753 if(variable == "pt") value = object->pt;
2754 else if(variable == "px") value = object->px;
2755 else if(variable == "py") value = object->py;
2756 else if(variable == "pz") value = object->pz;
2757 else if(variable == "phi") value = object->phi;
2758 else if(variable == "eta") value = object->eta;
2759 else if(variable == "theta") value = object->theta;
2760 else if(variable == "normChi2") value = object->normChi2;
2761 else if(variable == "dZ") value = object->dZ;
2762 else if(variable == "d0") value = object->d0;
2763 else if(variable == "d0err") value = object->d0err;
2764 else if(variable == "vx") value = object->vx;
2765 else if(variable == "vy") value = object->vy;
2766 else if(variable == "vz") value = object->vz;
2767 else if(variable == "charge") value = object->charge;
2768 else if(variable == "numValidHits") value = object->numValidHits;
2769 else if(variable == "isHighPurity") value = object->isHighPurity;
2770
2771 //additional BNs info for disappTrks
2772 else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
2773 else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
2774 else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
2775 else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
2776 else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
2777 else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
2778 else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
2779 else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
2780 else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
2781 else if(variable == "depTrkRp3") value = object->depTrkRp3;
2782 else if(variable == "trkRelIsoRp3") value = (object->depTrkRp3 - object->pt) / object->pt;
2783 else if(variable == "trkRelIsoRp5") value = (object->depTrkRp5 - object->pt) / object->pt;
2784 else if(variable == "depEcalRp3") value = object->depEcalRp3;
2785 else if(variable == "depHcalRp3") value = object->depHcalRp3;
2786 else if(variable == "depHoRp3") value = object->depHoRp3;
2787 else if(variable == "nTracksRp3") value = object->nTracksRp3;
2788 else if(variable == "trackerVetoPtRp3") value = object->trackerVetoPtRp3;
2789 else if(variable == "emVetoEtRp3") value = object->emVetoEtRp3;
2790 else if(variable == "hadVetoEtRp3") value = object->hadVetoEtRp3;
2791 else if(variable == "hoVetoEtRp3") value = object->hoVetoEtRp3;
2792 else if(variable == "depTrkRp5") value = object->depTrkRp5;
2793 else if(variable == "depEcalRp5") value = object->depEcalRp5;
2794 else if(variable == "depHcalRp5") value = object->depHcalRp5;
2795 else if(variable == "depHoRp5") value = object->depHoRp5;
2796 else if(variable == "nTracksRp5") value = object->nTracksRp5;
2797 else if(variable == "trackerVetoPtRp5") value = object->trackerVetoPtRp5;
2798 else if(variable == "emVetoEtRp5") value = object->emVetoEtRp5;
2799 else if(variable == "hadVetoEtRp5") value = object->hadVetoEtRp5;
2800 else if(variable == "hoVetoEtRp5") value = object->hoVetoEtRp5;
2801
2802 //user defined variables
2803 else if(variable == "d0wrtBS") value = (object->vx-events->at(0).BSx)*object->py/object->pt - (object->vy-events->at(0).BSy)*object->px/object->pt;
2804 else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
2805 else if(variable == "depTrkRp5MinusPt"){
2806 if ( (object->depTrkRp5 - object->pt) < 0 ) {
2807 // clog << "Warning: found track with depTrkRp5 < pt: depTrkRp5=" << object->depTrkRp5
2808 // << "; pt=" << object->pt
2809 // << "; object->depTrkRp5 - object->pt = " << object->depTrkRp5 - object->pt
2810 // << endl;
2811 value = 0;
2812 }
2813 else value = (object->depTrkRp5 - object->pt);
2814 }
2815 else if(variable == "depTrkRp3MinusPt"){
2816 if ( (object->depTrkRp3 - object->pt) < 0 ) {
2817 value = 0;
2818 }
2819 else value = (object->depTrkRp3 - object->pt);
2820 }
2821
2822 else if(variable == "dPhiMet") {
2823 if (const BNmet *met = chosenMET ()) {
2824 value = deltaPhi (object->phi, met->phi);
2825 } else value = -999;
2826 }
2827
2828
2829 else if(variable == "caloTotDeltaRp5") value = (object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
2830 else if(variable == "caloTotDeltaRp5ByP") value = ((object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
2831 else if(variable == "caloTotDeltaRp5RhoCorr") value = getTrkCaloTotRhoCorr(object);
2832 else if(variable == "caloTotDeltaRp5ByPRhoCorr") value = getTrkCaloTotRhoCorr(object) / pMag;
2833 else if(variable == "depTrkRp5RhoCorr") value = getTrkDepTrkRp5RhoCorr(object);
2834 else if(variable == "depTrkRp3RhoCorr") value = getTrkDepTrkRp3RhoCorr(object);
2835
2836 else if(variable == "depTrkRp5MinusPtRhoCorr") {
2837 if ( (getTrkDepTrkRp5RhoCorr(object) - object->pt) < 0 ) value = 0;
2838 else {value = (getTrkDepTrkRp5RhoCorr(object) - object->pt );}
2839 }
2840
2841 else if(variable == "depTrkRp3MinusPtRhoCorr")
2842 {
2843 if ( (getTrkDepTrkRp3RhoCorr(object) - object->pt) < 0 ) value = 0;
2844 else {value = (getTrkDepTrkRp3RhoCorr(object) - object->pt );}
2845 }
2846
2847 else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
2848 else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
2849 else if(variable == "ptErrorByPt") value = (object->ptError/object->pt);
2850 else if(variable == "ptError") value = object->ptError;
2851 else if(variable == "ptRes") value = getTrkPtRes(object);
2852 else if (variable == "d0wrtPV"){
2853 double vx = object->vx - chosenVertex ()->x,
2854 vy = object->vy - chosenVertex ()->y,
2855 px = object->px,
2856 py = object->py,
2857 pt = object->pt;
2858 value = (-vx * py + vy * px) / pt;
2859 }
2860 else if (variable == "dZwrtPV"){
2861 double vx = object->vx - chosenVertex ()->x,
2862 vy = object->vy - chosenVertex ()->y,
2863 vz = object->vz - chosenVertex ()->z,
2864 px = object->px,
2865 py = object->py,
2866 pz = object->pz,
2867 pt = object->pt;
2868 value = vz - (vx * px + vy * py)/pt * (pz/pt);
2869 }
2870
2871 else if(variable == "deltaRMinSubLeadJet") {
2872 // calculate minimum deltaR between track and any other subleading jet
2873 double trkJetDeltaRMin = 99.;
2874 for (uint ijet = 0; ijet<jets->size(); ijet++) {
2875 string empty = "";
2876 double isSubLeadingJet = valueLookup(&jets->at(ijet), "disappTrkSubLeadingJetID", "", empty);
2877 if (!isSubLeadingJet) continue; // only consider jets that pass the subleading jet ID criteria
2878 double jetEta = valueLookup(&jets->at(ijet), "eta", "", empty);
2879 double jetPhi = valueLookup(&jets->at(ijet), "phi", "", empty);
2880 double trkJetDeltaR = deltaR(object->eta, object->phi, jetEta, jetPhi);
2881 if (trkJetDeltaR < trkJetDeltaRMin) trkJetDeltaRMin = trkJetDeltaR;
2882 }
2883 value = trkJetDeltaRMin;
2884 }
2885
2886 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2887
2888 else if(variable == "genMatchedPdgId"){
2889 int index = getGenMatchedParticleIndex(object);
2890 if(index == -1) value = 0;
2891 else value = mcparticles->at(index).id;
2892 }
2893
2894
2895 else if(variable == "genMatchedId"){
2896 int index = getGenMatchedParticleIndex(object);
2897 if(index == -1) value = 0;
2898 else value = getPdgIdBinValue(mcparticles->at(index).id);
2899 }
2900 else if(variable == "genMatchedMotherId"){
2901 int index = getGenMatchedParticleIndex(object);
2902 if(index == -1) value = 0;
2903 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2904 }
2905 else if(variable == "genMatchedMotherIdReverse"){
2906 int index = getGenMatchedParticleIndex(object);
2907 if(index == -1) value = 24;
2908 else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2909 }
2910 else if(variable == "genMatchedGrandmotherId"){
2911 int index = getGenMatchedParticleIndex(object);
2912 if(index == -1) value = 0;
2913 else if(fabs(mcparticles->at(index).motherId) == 15){
2914 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2915 if(motherIndex == -1) value = 0;
2916 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2917 }
2918 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2919 }
2920 else if(variable == "genMatchedGrandmotherIdReverse"){
2921 int index = getGenMatchedParticleIndex(object);
2922 if(index == -1) value = 24;
2923 else if(fabs(mcparticles->at(index).motherId) == 15){
2924 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2925 if(motherIndex == -1) value = 24;
2926 else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2927 }
2928 else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2929 }
2930
2931
2932
2933 else{clog << "WARNING: invalid variable '" << variable << "' for BNtrack collection.\n"; value = -999;}
2934
2935 value = applyFunction(function, value);
2936
2937 return value;
2938 } // end track valueLookup
2939
2940
2941 //!genjet valueLookup
2942 double
2943 OSUAnalysis::valueLookup (const BNgenjet* object, string variable, string function, string &stringValue){
2944
2945 double value = 0.0;
2946
2947 if(variable == "pt") value = object->pt;
2948 else if(variable == "eta") value = object->eta;
2949 else if(variable == "phi") value = object->phi;
2950 else if(variable == "px") value = object->px;
2951 else if(variable == "py") value = object->py;
2952 else if(variable == "pz") value = object->pz;
2953 else if(variable == "et") value = object->et;
2954 else if(variable == "energy") value = object->energy;
2955 else if(variable == "mass") value = object->mass;
2956 else if(variable == "emEnergy") value = object->emEnergy;
2957 else if(variable == "hadEnergy") value = object->hadEnergy;
2958 else if(variable == "invisibleEnergy") value = object->invisibleEnergy;
2959 else if(variable == "auxiliaryEnergy") value = object->auxiliaryEnergy;
2960 else if(variable == "charge") value = object->charge;
2961
2962
2963 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2964
2965 value = applyFunction(function, value);
2966
2967 return value;
2968 }
2969
2970 //!mcparticle valueLookup
2971 double
2972 OSUAnalysis::valueLookup (const BNmcparticle* object, string variable, string function, string &stringValue){
2973
2974 double value = 0.0;
2975
2976 if(variable == "energy") value = object->energy;
2977 else if(variable == "et") value = object->et;
2978 else if(variable == "pt") value = object->pt;
2979 else if(variable == "px") value = object->px;
2980 else if(variable == "py") value = object->py;
2981 else if(variable == "pz") value = object->pz;
2982 else if(variable == "phi") value = object->phi;
2983 else if(variable == "eta") value = object->eta;
2984 else if(variable == "theta") value = object->theta;
2985 else if(variable == "mass") value = object->mass;
2986 else if(variable == "vx") value = object->vx;
2987 else if(variable == "vy") value = object->vy;
2988 else if(variable == "vz") value = object->vz;
2989 else if(variable == "motherET") value = object->motherET;
2990 else if(variable == "motherPT") value = object->motherPT;
2991 else if(variable == "motherPhi") value = object->motherPhi;
2992 else if(variable == "motherEta") value = object->motherEta;
2993 else if(variable == "mother0ET") value = object->mother0ET;
2994 else if(variable == "mother0PT") value = object->mother0PT;
2995 else if(variable == "mother0Phi") value = object->mother0Phi;
2996 else if(variable == "mother0Eta") value = object->mother0Eta;
2997 else if(variable == "mother1ET") value = object->mother1ET;
2998 else if(variable == "mother1PT") value = object->mother1PT;
2999 else if(variable == "mother1Phi") value = object->mother1Phi;
3000 else if(variable == "mother1Eta") value = object->mother1Eta;
3001 else if(variable == "daughter0ET") value = object->daughter0ET;
3002 else if(variable == "daughter0PT") value = object->daughter0PT;
3003 else if(variable == "daughter0Phi") value = object->daughter0Phi;
3004 else if(variable == "daughter0Eta") value = object->daughter0Eta;
3005 else if(variable == "daughter1ET") value = object->daughter1ET;
3006 else if(variable == "daughter1PT") value = object->daughter1PT;
3007 else if(variable == "daughter1Phi") value = object->daughter1Phi;
3008 else if(variable == "daughter1Eta") value = object->daughter1Eta;
3009 else if(variable == "grandMotherET") value = object->grandMotherET;
3010 else if(variable == "grandMotherPT") value = object->grandMotherPT;
3011 else if(variable == "grandMotherPhi") value = object->grandMotherPhi;
3012 else if(variable == "grandMotherEta") value = object->grandMotherEta;
3013 else if(variable == "grandMother00ET") value = object->grandMother00ET;
3014 else if(variable == "grandMother00PT") value = object->grandMother00PT;
3015 else if(variable == "grandMother00Phi") value = object->grandMother00Phi;
3016 else if(variable == "grandMother00Eta") value = object->grandMother00Eta;
3017 else if(variable == "grandMother01ET") value = object->grandMother01ET;
3018 else if(variable == "grandMother01PT") value = object->grandMother01PT;
3019 else if(variable == "grandMother01Phi") value = object->grandMother01Phi;
3020 else if(variable == "grandMother01Eta") value = object->grandMother01Eta;
3021 else if(variable == "grandMother10ET") value = object->grandMother10ET;
3022 else if(variable == "grandMother10PT") value = object->grandMother10PT;
3023 else if(variable == "grandMother10Phi") value = object->grandMother10Phi;
3024 else if(variable == "grandMother10Eta") value = object->grandMother10Eta;
3025 else if(variable == "grandMother11ET") value = object->grandMother11ET;
3026 else if(variable == "grandMother11PT") value = object->grandMother11PT;
3027 else if(variable == "grandMother11Phi") value = object->grandMother11Phi;
3028 else if(variable == "grandMother11Eta") value = object->grandMother11Eta;
3029 else if(variable == "charge") value = object->charge;
3030 else if(variable == "id") value = object->id;
3031 else if(variable == "status") value = object->status;
3032 else if(variable == "motherId") value = object->motherId;
3033 else if(variable == "motherCharge") value = object->motherCharge;
3034 else if(variable == "mother0Id") value = object->mother0Id;
3035 else if(variable == "mother0Status") value = object->mother0Status;
3036 else if(variable == "mother0Charge") value = object->mother0Charge;
3037 else if(variable == "mother1Id") value = object->mother1Id;
3038 else if(variable == "mother1Status") value = object->mother1Status;
3039 else if(variable == "mother1Charge") value = object->mother1Charge;
3040 else if(variable == "daughter0Id") value = object->daughter0Id;
3041 else if(variable == "daughter0Status") value = object->daughter0Status;
3042 else if(variable == "daughter0Charge") value = object->daughter0Charge;
3043 else if(variable == "daughter1Id") value = object->daughter1Id;
3044 else if(variable == "daughter1Status") value = object->daughter1Status;
3045 else if(variable == "daughter1Charge") value = object->daughter1Charge;
3046 else if(variable == "grandMotherId") value = object->grandMotherId;
3047 else if(variable == "grandMotherCharge") value = object->grandMotherCharge;
3048 else if(variable == "grandMother00Id") value = object->grandMother00Id;
3049 else if(variable == "grandMother00Status") value = object->grandMother00Status;
3050 else if(variable == "grandMother00Charge") value = object->grandMother00Charge;
3051 else if(variable == "grandMother01Id") value = object->grandMother01Id;
3052 else if(variable == "grandMother01Status") value = object->grandMother01Status;
3053 else if(variable == "grandMother01Charge") value = object->grandMother01Charge;
3054 else if(variable == "grandMother10Id") value = object->grandMother10Id;
3055 else if(variable == "grandMother10Status") value = object->grandMother10Status;
3056 else if(variable == "grandMother10Charge") value = object->grandMother10Charge;
3057 else if(variable == "grandMother11Id") value = object->grandMother11Id;
3058 else if(variable == "grandMother11Status") value = object->grandMother11Status;
3059 else if(variable == "grandMother11Charge") value = object->grandMother11Charge;
3060
3061 //user-defined variables
3062 else if (variable == "d0"){
3063 double vx = object->vx - chosenVertex ()->x,
3064 vy = object->vy - chosenVertex ()->y,
3065 px = object->px,
3066 py = object->py,
3067 pt = object->pt;
3068 value = (-vx * py + vy * px) / pt;
3069 }
3070 else if (variable == "dz"){
3071 double vx = object->vx - chosenVertex ()->x,
3072 vy = object->vy - chosenVertex ()->y,
3073 vz = object->vz - chosenVertex ()->z,
3074 px = object->px,
3075 py = object->py,
3076 pz = object->pz,
3077 pt = object->pt;
3078 value = vz - (vx * px + vy * py)/pt * (pz/pt);
3079 }
3080 else if(variable == "v0"){
3081 value = sqrt(object->vx*object->vx + object->vy*object->vy);
3082 }
3083 else if(variable == "deltaV0"){
3084 value = sqrt((object->vx-chosenVertex ()->x)*(object->vx-chosenVertex ()->x) + (object->vy-chosenVertex ()->y)*(object->vy-chosenVertex ()->y));
3085 }
3086 else if (variable == "deltaVx"){
3087 value = object->vx - chosenVertex ()->x;
3088 }
3089 else if (variable == "deltaVy"){
3090 value = object->vy - chosenVertex ()->y;
3091 }
3092 else if (variable == "deltaVz"){
3093 value = object->vz - chosenVertex ()->z;
3094 }
3095
3096
3097 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3098
3099 value = applyFunction(function, value);
3100
3101 return value;
3102 } // end mcparticle valueLookup
3103
3104
3105 //!primaryvertex valueLookup
3106 double
3107 OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable, string function, string &stringValue){
3108
3109 double value = 0.0;
3110
3111 if(variable == "x") value = object->x;
3112 else if(variable == "xError") value = object->xError;
3113 else if(variable == "y") value = object->y;
3114 else if(variable == "yError") value = object->yError;
3115 else if(variable == "z") value = object->z;
3116 else if(variable == "zError") value = object->zError;
3117 else if(variable == "rho") value = object->rho;
3118 else if(variable == "normalizedChi2") value = object->normalizedChi2;
3119 else if(variable == "ndof") value = object->ndof;
3120 else if(variable == "isFake") value = object->isFake;
3121 else if(variable == "isValid") value = object->isValid;
3122 else if(variable == "tracksSize") value = object->tracksSize;
3123 else if(variable == "isGood") value = object->isGood;
3124
3125
3126 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3127
3128 value = applyFunction(function, value);
3129
3130 return value;
3131 }
3132
3133 //!bxlumi valueLookup
3134 double
3135 OSUAnalysis::valueLookup (const BNbxlumi* object, string variable, string function, string &stringValue){
3136
3137 double value = 0.0;
3138
3139 if(variable == "bx_B1_now") value = object->bx_B1_now;
3140 else if(variable == "bx_B2_now") value = object->bx_B2_now;
3141 else if(variable == "bx_LUMI_now") value = object->bx_LUMI_now;
3142
3143
3144 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3145
3146 value = applyFunction(function, value);
3147
3148 return value;
3149 }
3150
3151 //!photon valueLookup
3152 double
3153 OSUAnalysis::valueLookup (const BNphoton* object, string variable, string function, string &stringValue){
3154
3155 double value = 0.0;
3156
3157 if(variable == "energy") value = object->energy;
3158 else if(variable == "et") value = object->et;
3159 else if(variable == "pt") value = object->pt;
3160 else if(variable == "px") value = object->px;
3161 else if(variable == "py") value = object->py;
3162 else if(variable == "pz") value = object->pz;
3163 else if(variable == "phi") value = object->phi;
3164 else if(variable == "eta") value = object->eta;
3165 else if(variable == "theta") value = object->theta;
3166 else if(variable == "trackIso") value = object->trackIso;
3167 else if(variable == "ecalIso") value = object->ecalIso;
3168 else if(variable == "hcalIso") value = object->hcalIso;
3169 else if(variable == "caloIso") value = object->caloIso;
3170 else if(variable == "trackIsoHollowConeDR03") value = object->trackIsoHollowConeDR03;
3171 else if(variable == "trackIsoSolidConeDR03") value = object->trackIsoSolidConeDR03;
3172 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
3173 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
3174 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
3175 else if(variable == "trackIsoHollowConeDR04") value = object->trackIsoHollowConeDR04;
3176 else if(variable == "trackIsoSolidConeDR04") value = object->trackIsoSolidConeDR04;
3177 else if(variable == "ecalIsoDR04") value = object->ecalIsoDR04;
3178 else if(variable == "hcalIsoDR04") value = object->hcalIsoDR04;
3179 else if(variable == "caloIsoDR04") value = object->caloIsoDR04;
3180 else if(variable == "hadOverEm") value = object->hadOverEm;
3181 else if(variable == "sigmaEtaEta") value = object->sigmaEtaEta;
3182 else if(variable == "sigmaIetaIeta") value = object->sigmaIetaIeta;
3183 else if(variable == "r9") value = object->r9;
3184 else if(variable == "scEnergy") value = object->scEnergy;
3185 else if(variable == "scRawEnergy") value = object->scRawEnergy;
3186 else if(variable == "scSeedEnergy") value = object->scSeedEnergy;
3187 else if(variable == "scEta") value = object->scEta;
3188 else if(variable == "scPhi") value = object->scPhi;
3189 else if(variable == "scZ") value = object->scZ;
3190 else if(variable == "genET") value = object->genET;
3191 else if(variable == "genPT") value = object->genPT;
3192 else if(variable == "genPhi") value = object->genPhi;
3193 else if(variable == "genEta") value = object->genEta;
3194 else if(variable == "genMotherET") value = object->genMotherET;
3195 else if(variable == "genMotherPT") value = object->genMotherPT;
3196 else if(variable == "genMotherPhi") value = object->genMotherPhi;
3197 else if(variable == "genMotherEta") value = object->genMotherEta;
3198 else if(variable == "eMax") value = object->eMax;
3199 else if(variable == "eLeft") value = object->eLeft;
3200 else if(variable == "eRight") value = object->eRight;
3201 else if(variable == "eTop") value = object->eTop;
3202 else if(variable == "eBottom") value = object->eBottom;
3203 else if(variable == "e3x3") value = object->e3x3;
3204 else if(variable == "swissCross") value = object->swissCross;
3205 else if(variable == "seedEnergy") value = object->seedEnergy;
3206 else if(variable == "seedTime") value = object->seedTime;
3207 else if(variable == "swissCrossNoI85") value = object->swissCrossNoI85;
3208 else if(variable == "swissCrossI85") value = object->swissCrossI85;
3209 else if(variable == "E2overE9NoI85") value = object->E2overE9NoI85;
3210 else if(variable == "E2overE9I85") value = object->E2overE9I85;
3211 else if(variable == "IDTight") value = object->IDTight;
3212 else if(variable == "IDLoose") value = object->IDLoose;
3213 else if(variable == "IDLooseEM") value = object->IDLooseEM;
3214 else if(variable == "genId") value = object->genId;
3215 else if(variable == "genCharge") value = object->genCharge;
3216 else if(variable == "genMotherId") value = object->genMotherId;
3217 else if(variable == "genMotherCharge") value = object->genMotherCharge;
3218 else if(variable == "isEB") value = object->isEB;
3219 else if(variable == "isEE") value = object->isEE;
3220 else if(variable == "isGap") value = object->isGap;
3221 else if(variable == "isEBEEGap") value = object->isEBEEGap;
3222 else if(variable == "isEBGap") value = object->isEBGap;
3223 else if(variable == "isEEGap") value = object->isEEGap;
3224 else if(variable == "hasPixelSeed") value = object->hasPixelSeed;
3225 else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
3226
3227
3228
3229
3230 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
3231
3232 else if(variable == "genMatchedPdgId"){
3233 int index = getGenMatchedParticleIndex(object);
3234 if(index == -1) value = 0;
3235 else value = mcparticles->at(index).id;
3236 }
3237
3238
3239
3240 else if(variable == "genMatchedId"){
3241 int index = getGenMatchedParticleIndex(object);
3242 if(index == -1) value = 0;
3243 else value = getPdgIdBinValue(mcparticles->at(index).id);
3244 }
3245 else if(variable == "genMatchedMotherId"){
3246 int index = getGenMatchedParticleIndex(object);
3247 if(index == -1) value = 0;
3248 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
3249 }
3250 else if(variable == "genMatchedMotherIdReverse"){
3251 int index = getGenMatchedParticleIndex(object);
3252 if(index == -1) value = 24;
3253 else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
3254 }
3255 else if(variable == "genMatchedGrandmotherId"){
3256 int index = getGenMatchedParticleIndex(object);
3257 if(index == -1) value = 0;
3258 else if(fabs(mcparticles->at(index).motherId) == 15){
3259 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
3260 if(motherIndex == -1) value = 0;
3261 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
3262 }
3263 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
3264 }
3265 else if(variable == "genMatchedGrandmotherIdReverse"){
3266 int index = getGenMatchedParticleIndex(object);
3267 if(index == -1) value = 24;
3268 else if(fabs(mcparticles->at(index).motherId) == 15){
3269 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
3270 if(motherIndex == -1) value = 24;
3271 else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
3272 }
3273 else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
3274 }
3275
3276
3277 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3278
3279 value = applyFunction(function, value);
3280
3281 return value;
3282 } // end photon valueLookup
3283
3284
3285 //!supercluster valueLookup
3286 double
3287 OSUAnalysis::valueLookup (const BNsupercluster* object, string variable, string function, string &stringValue){
3288
3289 double value = 0.0;
3290
3291 if(variable == "energy") value = object->energy;
3292 else if(variable == "et") value = object->et;
3293 else if(variable == "ex") value = object->ex;
3294 else if(variable == "ey") value = object->ey;
3295 else if(variable == "ez") value = object->ez;
3296 else if(variable == "phi") value = object->phi;
3297 else if(variable == "eta") value = object->eta;
3298 else if(variable == "theta") value = object->theta;
3299
3300
3301 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3302
3303 value = applyFunction(function, value);
3304
3305 return value;
3306 }
3307
3308 //!trigobj valueLookup
3309 double
3310 OSUAnalysis::valueLookup (const BNtrigobj* object, string variable, string function, string &stringValue){
3311
3312 double value = 0.0;
3313
3314 if(variable == "pt") value = object->pt;
3315 else if(variable == "eta") value = object->eta;
3316 else if(variable == "phi") value = object->phi;
3317 else if(variable == "px") value = object->px;
3318 else if(variable == "py") value = object->py;
3319 else if(variable == "pz") value = object->pz;
3320 else if(variable == "et") value = object->et;
3321 else if(variable == "energy") value = object->energy;
3322 else if(variable == "etTotal") value = object->etTotal;
3323 else if(variable == "id") value = object->id;
3324 else if(variable == "charge") value = object->charge;
3325 else if(variable == "isIsolated") value = object->isIsolated;
3326 else if(variable == "isMip") value = object->isMip;
3327 else if(variable == "isForward") value = object->isForward;
3328 else if(variable == "isRPC") value = object->isRPC;
3329 else if(variable == "bx") value = object->bx;
3330 else if(variable == "filter") {
3331 if ((stringValue = object->filter) == "")
3332 stringValue = "none"; // stringValue should only be empty if value is filled
3333 }
3334
3335 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3336
3337 value = applyFunction(function, value);
3338
3339 return value;
3340 }
3341
3342 //!muon-muon pair valueLookup
3343 double
3344 OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function, string &stringValue){
3345
3346 double value = 0.0;
3347
3348 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3349 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3350 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3351 else if(variable == "invMass"){
3352 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3353 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3354 value = (fourVector1 + fourVector2).M();
3355 }
3356 else if(variable == "pt"){
3357 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3358 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3359 value = (fourVector1 + fourVector2).Pt();
3360 }
3361 else if(variable == "threeDAngle")
3362 {
3363 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3364 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3365 value = (threeVector1.Angle(threeVector2));
3366 }
3367 else if(variable == "alpha")
3368 {
3369 static const double pi = 3.1415926535897932384626433832795028841971693993751058;
3370 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3371 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3372 value = (pi-threeVector1.Angle(threeVector2));
3373 }
3374 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
3375 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
3376 else if(variable == "d0Sign"){
3377 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
3378 if(d0Sign < 0) value = -0.5;
3379 else if (d0Sign > 0) value = 0.5;
3380 else value = -999;
3381 }
3382 else if(variable == "chargeProduct"){
3383 value = object1->charge*object2->charge;
3384 }
3385 else if(variable == "muon1CorrectedD0Vertex"){
3386 value = object1->correctedD0Vertex;
3387 }
3388 else if(variable == "muon2CorrectedD0Vertex"){
3389 value = object2->correctedD0Vertex;
3390 }
3391 else if(variable == "muon1timeAtIpInOut"){
3392 value = object1->timeAtIpInOut;
3393 }
3394 else if(variable == "muon2timeAtIpInOut"){
3395 value = object2->timeAtIpInOut;
3396 }
3397 else if(variable == "muon1correctedD0")
3398 {
3399 value = object1->correctedD0;
3400 }
3401 else if(variable == "muon2correctedD0")
3402 {
3403 value = object2->correctedD0;
3404 }
3405
3406 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3407
3408 value = applyFunction(function, value);
3409
3410 return value;
3411 } // end muon-muon pair valueLookup
3412
3413
3414 //!muon-photon pair valueLookup
3415 double
3416 OSUAnalysis::valueLookup (const BNmuon* object1, const BNphoton* object2, string variable, string function, string &stringValue){
3417
3418 double value = 0.0;
3419
3420 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3421 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3422 else if(variable == "photonEta") value = object2->eta;
3423 else if(variable == "photonPt") value = object2->pt;
3424 else if(variable == "muonEta") value = object1->eta;
3425 else if(variable == "photonPhi") value = object2->phi;
3426 else if(variable == "muonPhi") value = object1->phi;
3427 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3428 else if(variable == "photonGenMotherId") value = object2->genMotherId;
3429 else if(variable == "muonRelPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3430 else if(variable == "invMass"){
3431 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3432 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3433 value = (fourVector1 + fourVector2).M();
3434 }
3435 else if(variable == "pt"){
3436 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3437 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3438 value = (fourVector1 + fourVector2).Pt();
3439 }
3440 else if(variable == "threeDAngle")
3441 {
3442 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3443 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3444 value = (threeVector1.Angle(threeVector2));
3445 }
3446 else if(variable == "alpha")
3447 {
3448 static const double pi = 3.1415926535897932384626433832795028841971693993751058;
3449 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3450 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3451 value = (pi-threeVector1.Angle(threeVector2));
3452 }
3453 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3454
3455 value = applyFunction(function, value);
3456
3457 return value;
3458 }
3459
3460 //!electron-photon pair valueLookup
3461 double
3462 OSUAnalysis::valueLookup (const BNelectron* object1, const BNphoton* object2, string variable, string function, string &stringValue){
3463
3464 double value = 0.0;
3465
3466 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3467 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3468 else if(variable == "photonEta") value = object2->eta;
3469 else if(variable == "photonPt") value = object2->pt;
3470 else if(variable == "electronEta") value = object1->eta;
3471 else if(variable == "photonPhi") value = object2->phi;
3472 else if(variable == "electronPhi") value = object1->phi;
3473 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3474 else if(variable == "photonGenMotherId") value = object2->genMotherId;
3475 else if(variable == "electronRelPFrhoIso") value = ( object1->chargedHadronIsoDR03 + max(0.0, object1->neutralHadronIsoDR03 + object1->photonIsoDR03 - object1->AEffDr03*object1->rhoPrime) ) / object1->pt;
3476 else if(variable == "invMass"){
3477 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3478 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3479 value = (fourVector1 + fourVector2).M();
3480 }
3481 else if(variable == "pt"){
3482 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3483 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3484 value = (fourVector1 + fourVector2).Pt();
3485 }
3486 else if(variable == "threeDAngle")
3487 {
3488 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3489 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3490 value = (threeVector1.Angle(threeVector2));
3491 }
3492 else if(variable == "alpha")
3493 {
3494 static const double pi = 3.1415926535897932384626433832795028841971693993751058;
3495 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3496 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3497 value = (pi-threeVector1.Angle(threeVector2));
3498 }
3499 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3500
3501 value = applyFunction(function, value);
3502
3503 return value;
3504 }
3505
3506 //!electron-electron pair valueLookup
3507 double
3508 OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function, string &stringValue){
3509
3510 double value = 0.0;
3511
3512 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3513 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3514 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3515 else if(variable == "invMass"){
3516 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3517 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3518 value = (fourVector1 + fourVector2).M();
3519 }
3520 else if(variable == "pt"){
3521 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3522 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3523 value = (fourVector1 + fourVector2).Pt();
3524 }
3525 else if(variable == "threeDAngle")
3526 {
3527 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3528 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3529 value = (threeVector1.Angle(threeVector2));
3530 }
3531 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
3532 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
3533 else if(variable == "d0Sign"){
3534 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
3535 if(d0Sign < 0) value = -0.5;
3536 else if (d0Sign > 0) value = 0.5;
3537 else value = -999;
3538 }
3539 else if(variable == "chargeProduct"){
3540 value = object1->charge*object2->charge;
3541 }
3542 else if(variable == "electron1CorrectedD0Vertex"){
3543 value = object1->correctedD0Vertex;
3544 }
3545 else if(variable == "electron2CorrectedD0Vertex"){
3546 value = object2->correctedD0Vertex;
3547 }
3548 else if(variable == "electron1CorrectedD0"){
3549 value = object1->correctedD0;
3550 }
3551 else if(variable == "electron2CorrectedD0"){
3552 value = object2->correctedD0;
3553 }
3554
3555 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3556
3557 value = applyFunction(function, value);
3558
3559 return value;
3560 }
3561
3562 //!electron-muon pair valueLookup
3563 double
3564 OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function, string &stringValue){
3565
3566 double value = 0.0;
3567
3568 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3569 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3570 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3571 else if(variable == "invMass"){
3572 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3573 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3574 value = (fourVector1 + fourVector2).M();
3575 }
3576 else if(variable == "pt"){
3577 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3578 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3579 value = (fourVector1 + fourVector2).Pt();
3580 }
3581 else if(variable == "threeDAngle")
3582 {
3583 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3584 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3585 value = (threeVector1.Angle(threeVector2));
3586 }
3587 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
3588 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
3589 else if(variable == "d0Sign"){
3590 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
3591 if(d0Sign < 0) value = -0.5;
3592 else if (d0Sign > 0) value = 0.5;
3593 else value = -999;
3594 }
3595 else if(variable == "chargeProduct"){
3596 value = object1->charge*object2->charge;
3597 }
3598 else if(variable == "electronCorrectedD0Vertex"){
3599 value = object1->correctedD0Vertex;
3600 }
3601 else if(variable == "muonCorrectedD0Vertex"){
3602 value = object2->correctedD0Vertex;
3603 }
3604 else if(variable == "electronCorrectedD0"){
3605 value = object1->correctedD0;
3606 }
3607 else if(variable == "muonCorrectedD0"){
3608 value = object2->correctedD0;
3609 }
3610 else if(variable == "electronDetIso"){
3611 value = (object1->trackIso) / object1->pt;
3612 }
3613 else if(variable == "muonDetIso"){
3614 value = (object2->trackIsoDR03) / object2->pt;
3615 }
3616 else if(variable == "electronRelPFrhoIso"){
3617 value = ( object1->chargedHadronIsoDR03 + max(0.0, object1->neutralHadronIsoDR03 + object1->photonIsoDR03 - object1->AEffDr03*object1->rhoPrime) ) / object1->pt;
3618 }
3619 else if(variable == "muonRelPFdBetaIso"){
3620 value = (object2->pfIsoR04SumChargedHadronPt + max(0.0, object2->pfIsoR04SumNeutralHadronEt + object2->pfIsoR04SumPhotonEt - 0.5*object2->pfIsoR04SumPUPt)) / object2->pt;
3621 }
3622 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3623 value = applyFunction(function, value);
3624
3625 return value;
3626 } // end electron-muon pair valueLookup
3627
3628
3629 //!electron-jet pair valueLookup
3630 double
3631 OSUAnalysis::valueLookup (const BNelectron* object1, const BNjet* object2, string variable, string function, string &stringValue){
3632
3633 double value = 0.0;
3634
3635 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3636 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3637 else if(variable == "jetEta") value = object2->eta;
3638 else if(variable == "jetPhi") value = object2->phi;
3639 else if(variable == "electronEta") value = object1->eta;
3640 else if(variable == "electronPhi") value = object1->phi;
3641 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3642 else if(variable == "invMass"){
3643 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3644 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3645 value = (fourVector1 + fourVector2).M();
3646 }
3647 else if(variable == "pt"){
3648 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3649 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3650 value = (fourVector1 + fourVector2).Pt();
3651 }
3652 else if(variable == "threeDAngle")
3653 {
3654 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3655 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3656 value = (threeVector1.Angle(threeVector2));
3657 }
3658 else if(variable == "chargeProduct"){
3659 value = object1->charge*object2->charge;
3660 }
3661
3662 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3663 value = applyFunction(function, value);
3664
3665 return value;
3666 }
3667
3668 //!photon-jet pair valueLookup
3669 double
3670 OSUAnalysis::valueLookup (const BNphoton* object1, const BNjet* object2, string variable, string function, string &stringValue){
3671
3672 double value = 0.0;
3673
3674 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3675 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3676 else if(variable == "jetEta") value = object2->eta;
3677 else if(variable == "jetPhi") value = object2->phi;
3678 else if(variable == "photonEta") value = object1->eta;
3679 else if(variable == "photonPhi") value = object1->phi;
3680 else if(variable == "photonGenMotherId") value = object1->genMotherId;
3681 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3682 else if(variable == "invMass"){
3683 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3684 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3685 value = (fourVector1 + fourVector2).M();
3686 }
3687 else if(variable == "pt"){
3688 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3689 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3690 value = (fourVector1 + fourVector2).Pt();
3691 }
3692 else if(variable == "threeDAngle")
3693 {
3694 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3695 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3696 value = (threeVector1.Angle(threeVector2));
3697 }
3698 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3699 value = applyFunction(function, value);
3700
3701 return value;
3702 }
3703
3704 // track-jet pair valueLookup
3705 double
3706 OSUAnalysis::valueLookup (const BNtrack* object1, const BNjet* object2, string variable, string function, string &stringValue){
3707
3708 double value = 0.0;
3709
3710 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3711 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3712
3713 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3714 value = applyFunction(function, value);
3715
3716 return value;
3717
3718 }
3719
3720
3721
3722 // met-jet pair valueLookup
3723 double
3724 OSUAnalysis::valueLookup (const BNmet* object1, const BNjet* object2, string variable, string function, string &stringValue){
3725
3726 double value = 0.0;
3727
3728 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3729
3730 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3731 value = applyFunction(function, value);
3732
3733 return value;
3734
3735 }
3736
3737
3738
3739 //!muon-jet pair valueLookup
3740 double
3741 OSUAnalysis::valueLookup (const BNmuon* object1, const BNjet* object2, string variable, string function, string &stringValue){
3742
3743 double value = 0.0;
3744
3745 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3746 else if(variable == "jetEta") value = object2->eta;
3747 else if(variable == "relPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3748 else if(variable == "jetPt") value = object2->pt;
3749 else if(variable == "jetPhi") value = object2->phi;
3750 else if(variable == "deltaPt") value = object1->pt - object2->pt;
3751 else if(variable == "muonEta") value = object1->eta;
3752 else if(variable == "muonPt") value = object1->pt;
3753 else if(variable == "muonPhi") value = object1->phi;
3754 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3755 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3756 else if(variable == "invMass"){
3757 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3758 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3759 value = (fourVector1 + fourVector2).M();
3760 }
3761 else if(variable == "pt"){
3762 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3763 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3764 value = (fourVector1 + fourVector2).Pt();
3765 }
3766 else if(variable == "threeDAngle")
3767 {
3768 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3769 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3770 value = (threeVector1.Angle(threeVector2));
3771 }
3772 else if(variable == "chargeProduct"){
3773 value = object1->charge*object2->charge;
3774 }
3775
3776 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3777 value = applyFunction(function, value);
3778
3779 return value;
3780 }
3781
3782 //!muon-event valueLookup
3783 double
3784 OSUAnalysis::valueLookup (const BNmuon* object1, const BNevent* object2, string variable, string function, string &stringValue){
3785
3786 double value = 0.0;
3787
3788 if(variable == "muonEta") value = object1->eta;
3789 else if(variable == "muonPt") value = object1->pt;
3790 else if(variable == "muonPhi") value = object1->phi;
3791 else if(variable == "Ht") value = getHt(jets.product());
3792 else if(variable == "pthat") value = object2->pthat;
3793 else if(variable == "relPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3794 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3795 value = applyFunction(function, value);
3796
3797 return value;
3798 }
3799 //!jet-jet pair valueLookup
3800 double
3801 OSUAnalysis::valueLookup (const BNjet* object1, const BNjet* object2, string variable, string function, string &stringValue){
3802
3803 double value = 0.0;
3804
3805 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3806 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3807 else if(variable == "deltaPt") value = object1->pt - object2->pt;
3808 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3809 else if(variable == "invMass"){
3810 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3811 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3812 value = (fourVector1 + fourVector2).M();
3813 }
3814 else if(variable == "pt"){
3815 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3816 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3817 value = (fourVector1 + fourVector2).Pt();
3818 }
3819 else if(variable == "threeDAngle")
3820 {
3821 TVector3 threeVector1(object1->px, object1->py, object1->pz);
3822 TVector3 threeVector2(object2->px, object2->py, object2->pz);
3823 value = (threeVector1.Angle(threeVector2));
3824 }
3825 else if(variable == "chargeProduct"){
3826 value = object1->charge*object2->charge;
3827 }
3828
3829 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3830 value = applyFunction(function, value);
3831
3832 return value;
3833 }
3834
3835 //!electron-track pair valueLookup
3836 double
3837 OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3838 double electronMass = 0.000511;
3839 double value = 0.0;
3840 TLorentzVector fourVector1(0, 0, 0, 0);
3841 TLorentzVector fourVector2(0, 0, 0, 0);
3842 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3843 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3844 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3845 else if(variable == "invMass"){
3846 fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, electronMass);
3847 fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, electronMass );
3848
3849 value = (fourVector1 + fourVector2).M();
3850 }
3851 else if(variable == "chargeProduct"){
3852 value = object1->charge*object2->charge;
3853 }
3854 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3855 value = applyFunction(function, value);
3856 return value;
3857
3858 }
3859
3860
3861 //!muon-track pair valueLookup
3862 double
3863 OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3864 double pionMass = 0.140;
3865 double muonMass = 0.106;
3866 double value = 0.0;
3867 TLorentzVector fourVector1(0, 0, 0, 0);
3868 TLorentzVector fourVector2(0, 0, 0, 0);
3869 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3870 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3871 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3872 else if(variable == "invMass"){
3873 fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, muonMass);
3874 fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, pionMass );
3875
3876 value = (fourVector1 + fourVector2).M();
3877 }
3878 else if(variable == "chargeProduct"){
3879 value = object1->charge*object2->charge;
3880 }
3881
3882 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3883 value = applyFunction(function, value);
3884 return value;
3885 }
3886
3887 //!tau-tau pair valueLookup
3888 double
3889 OSUAnalysis::valueLookup (const BNtau* object1, const BNtau* object2, string variable, string function, string &stringValue){
3890 double value = 0.0;
3891 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3892 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3893 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3894 else if(variable == "invMass"){
3895 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3896 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3897 value = (fourVector1 + fourVector2).M();
3898 }
3899
3900 else if(variable == "chargeProduct"){
3901 value = object1->charge*object2->charge;
3902 }
3903
3904 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3905 value = applyFunction(function, value);
3906 return value;
3907 }
3908
3909 //!muon-tau pair valueLookup
3910 double
3911 OSUAnalysis::valueLookup (const BNmuon* object1, const BNtau* object2, string variable, string function, string &stringValue){
3912 double value = 0.0;
3913 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3914 else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3915 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3916 else if(variable == "invMass"){
3917 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3918 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3919 value = (fourVector1 + fourVector2).M();
3920 }
3921
3922 else if(variable == "chargeProduct"){
3923 value = object1->charge*object2->charge;
3924 }
3925
3926 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3927 value = applyFunction(function, value);
3928 return value;
3929 }
3930
3931 //!tau-track pair valueLookup
3932 double
3933 OSUAnalysis::valueLookup (const BNtau* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3934 double value = 0.0;
3935 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3936 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3937 else if(variable == "chargeProduct"){
3938 value = object1->charge*object2->charge;
3939 }
3940
3941 else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3942 value = applyFunction(function, value);
3943 return value;
3944 }
3945
3946
3947 //!track-event pair valueLookup
3948 double
3949 OSUAnalysis::valueLookup (const BNtrack* object1, const BNevent* object2, string variable, string function, string &stringValue){
3950
3951 double value = 0.0;
3952 double pMag = sqrt(object1->pt * object1->pt +
3953 object1->pz * object1->pz);
3954
3955 if (variable == "numPV") value = object2->numPV;
3956 else if (variable == "caloTotDeltaRp5") value = (object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5);
3957 else if (variable == "caloTotDeltaRp5ByP") value = ((object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5)/pMag);
3958 else if (variable == "caloTotDeltaRp5_RhoCorr") value = getTrkCaloTotRhoCorr(object1);
3959 else if (variable == "caloTotDeltaRp5ByP_RhoCorr") value = getTrkCaloTotRhoCorr(object1) / pMag;
3960
3961 else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3962
3963 value = applyFunction(function, value);
3964
3965 return value;
3966
3967 }
3968
3969 //!electron-trigobj pair valueLookup
3970 double
3971 OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrigobj* object2, string variable, string function, string &stringValue){
3972
3973 double value = 0.0;
3974
3975 if (variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3976 else if (variable == "match"){
3977 if (deltaR(object1->eta,object1->phi,object2->eta,object2->phi) < 0.2 && abs(object2->id) == 11)
3978 stringValue = object2->filter;
3979 else
3980 stringValue = "none";
3981 }
3982
3983 else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3984
3985 value = applyFunction(function, value);
3986
3987 return value;
3988
3989 }
3990
3991 //!muon-trigobj pair valueLookup
3992 double
3993 OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrigobj* object2, string variable, string function, string &stringValue){
3994
3995 double value = 0.0;
3996
3997 if (variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3998
3999 else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
4000
4001 value = applyFunction(function, value);
4002
4003 return value;
4004
4005 }
4006
4007 //!stop valueLookup
4008 double
4009 OSUAnalysis::valueLookup (const BNstop* object, string variable, string function, string &stringValue){
4010
4011
4012 double value = 0.0;
4013
4014 if(variable == "ctau") value = object->ctau;
4015
4016 else if (variable == "d0"){
4017 double vx = object->vx - chosenVertex ()->x,
4018 vy = object->vy - chosenVertex ()->y,
4019 px = object->px,
4020 py = object->py,
4021 pt = object->pt;
4022 value = (-vx * py + vy * px) / pt;
4023 }
4024
4025 else if (variable == "dz"){
4026 double vx = object->vx - chosenVertex ()->x,
4027 vy = object->vy - chosenVertex ()->y,
4028 vz = object->vz - chosenVertex ()->z,
4029 px = object->px,
4030 py = object->py,
4031 pz = object->pz,
4032 pt = object->pt;
4033 value = vz - (vx * px + vy * py)/pt * (pz/pt);
4034 }
4035
4036 else if (variable == "minD0"){
4037 double minD0=999;
4038 for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4039 double vx = object->vx - vertex->x,
4040 vy = object->vy - vertex->y,
4041 px = object->px,
4042 py = object->py,
4043 pt = object->pt;
4044 value = (-vx * py + vy * px) / pt;
4045 if(abs(value) < abs(minD0)) minD0 = value;
4046 }
4047 value = minD0;
4048 }
4049 else if (variable == "minDz"){
4050 double minDz=999;
4051 for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4052 double vx = object->vx - vertex->x,
4053 vy = object->vy - vertex->y,
4054 vz = object->vz - vertex->z,
4055 px = object->px,
4056 py = object->py,
4057 pz = object->pz,
4058 pt = object->pt;
4059 value = vz - (vx * px + vy * py)/pt * (pz/pt);
4060 if(abs(value) < abs(minDz)) minDz = value;
4061 }
4062 value = minDz;
4063 }
4064 else if(variable == "distToVertex"){
4065 value = sqrt((object->vx-chosenVertex()->x)*(object->vx-chosenVertex()->x) + \
4066 (object->vy-chosenVertex()->y)*(object->vy-chosenVertex()->y) + \
4067 (object->vz-chosenVertex()->z)*(object->vz-chosenVertex()->z));
4068 }
4069 else if (variable == "minDistToVertex"){
4070 double minDistToVertex=999;
4071 for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4072 value = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4073 (object->vy-vertex->y)*(object->vy-vertex->y) + \
4074 (object->vz-vertex->z)*(object->vz-vertex->z));
4075
4076 if(abs(value) < abs(minDistToVertex)) minDistToVertex = value;
4077 }
4078 value = minDistToVertex;
4079 }
4080 else if (variable == "distToVertexDifference"){
4081 double minDistToVertex=999;
4082 for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4083 value = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4084 (object->vy-vertex->y)*(object->vy-vertex->y) + \
4085 (object->vz-vertex->z)*(object->vz-vertex->z));
4086
4087 if(abs(value) < abs(minDistToVertex)) minDistToVertex = value;
4088 }
4089 double distToChosenVertex = sqrt((object->vx-chosenVertex()->x)*(object->vx-chosenVertex()->x) + \
4090 (object->vy-chosenVertex()->y)*(object->vy-chosenVertex()->y) + \
4091 (object->vz-chosenVertex()->z)*(object->vz-chosenVertex()->z));
4092
4093 value = distToChosenVertex - minDistToVertex;
4094 }
4095
4096 else if (variable == "closestVertexRank"){
4097 double minDistToVertex=999;
4098 int vertex_rank = 0;
4099 for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4100 vertex_rank++;
4101 int dist = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4102 (object->vy-vertex->y)*(object->vy-vertex->y) + \
4103 (object->vz-vertex->z)*(object->vz-vertex->z));
4104
4105 if(abs(dist) < abs(minDistToVertex)){
4106 value = vertex_rank;
4107 minDistToVertex = dist;
4108 }
4109 }
4110 }
4111
4112 else if (variable == "decaysToTau"){
4113 value = abs (object->daughter0Id) == 15 || abs (object->daughter1Id) == 15;
4114 }
4115
4116
4117
4118
4119 else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
4120
4121 value = applyFunction(function, value);
4122
4123 return value;
4124
4125 } // end stop valueLookup
4126
4127
4128
4129
4130
4131 // Calculate the number of tracks in cone of DeltaR<0.5 around track1.
4132 // Return true iff no other tracks are found in this cone.
4133 int
4134 OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
4135 for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
4136 if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.
4137 double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
4138 if(deltaRtrk < 0.5) return 0;
4139 }
4140 return 1;
4141
4142 }
4143
4144 //calculate the scalar sum of Jet Pt in the event.
4145 double
4146 OSUAnalysis::getHt (const BNjetCollection* jetColl){
4147 double Ht = 0;
4148 for(BNjetCollection::const_iterator jet = jetColl->begin(); jet !=jetColl->end(); jet++){
4149 Ht += abs(jet->pt);
4150 }
4151 return Ht;
4152 }
4153
4154 double
4155 OSUAnalysis::getTrkPtRes (const BNtrack* track1){
4156
4157 double ptTrue = getTrkPtTrue(track1, mcparticles.product());
4158 double PtRes = (track1->pt - ptTrue) / ptTrue;
4159
4160 return PtRes;
4161
4162 }
4163
4164
4165 double
4166 OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
4167 double value = -99;
4168 double genDeltaRLowest = 999;
4169
4170 for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
4171 double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
4172 if (genDeltaRtemp < genDeltaRLowest) {
4173 genDeltaRLowest = genDeltaRtemp;
4174 if (genDeltaRLowest < 0.05) { // Only consider it truth-matched if DeltaR<0.15.
4175 double ptTrue = genPart->pt;
4176 value = ptTrue;
4177 }
4178 }
4179 }
4180
4181 return value;
4182
4183 }
4184
4185 double
4186 OSUAnalysis::getTrkCaloTotRhoCorr(const BNtrack* track) {
4187 // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.
4188 if (!useTrackCaloRhoCorr_) return -99;
4189 // if (!rhokt6CaloJetsHandle_) {
4190 // clog << "ERROR [getTrkCaloTotRhoCorr]: The collection rhokt6CaloJetsHandle is not available!" << endl;
4191 // return -99;
4192 // }
4193 double radDeltaRCone = 0.5;
4194 double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2); // Define effective area as pi*r^2, where r is radius of DeltaR cone.
4195 double rawCaloTot = track->caloHadDeltaRp5 + track->caloEMDeltaRp5;
4196 double caloTotRhoCorrCalo = TMath::Max(0., rawCaloTot - rhoCorr_kt6CaloJets);
4197 return caloTotRhoCorrCalo;
4198
4199 }
4200
4201 double
4202 OSUAnalysis::getTrkDepTrkRp5RhoCorr(const BNtrack* track) {
4203 // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.
4204 if (!useTrackCaloRhoCorr_) return -99;
4205 // if (!rhokt6CaloJetsHandle_) {
4206 // clog << "ERROR [getTrkCaloTotRhoCorr]: The collection rhokt6CaloJetsHandle is not available!" << endl;
4207 // return -99;
4208 // }
4209 double radDeltaRCone = 0.5;
4210 double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2); // Define effective area as pi*r^2, where r is radius of DeltaR cone.
4211 double rawDepTrkRp5 = track->depTrkRp5;
4212 double depTrkRp5RhoCorr = TMath::Max(0., rawDepTrkRp5 - rhoCorr_kt6CaloJets);
4213 return depTrkRp5RhoCorr;
4214
4215 }
4216
4217 double
4218 OSUAnalysis::getTrkDepTrkRp3RhoCorr(const BNtrack* track) {
4219 // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.
4220 if (!useTrackCaloRhoCorr_) return -99;
4221 // if (!rhokt6CaloJetsHandle_) {
4222 // clog << "ERROR [getTrkCaloTotRhoCorr]: The collection rhokt6CaloJetsHandle is not available!" << endl;
4223 // return -99;
4224 // }
4225 double radDeltaRCone = 0.3;
4226 // Define effective area as pi*r^2, where r is radius of DeltaR cone
4227 double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);
4228 double rawDepTrkRp3 = track->depTrkRp3;
4229 double depTrkRp3RhoCorr = TMath::Max(0., rawDepTrkRp3 - rhoCorr_kt6CaloJets);
4230 return depTrkRp3RhoCorr;
4231
4232 }
4233
4234
4235
4236
4237 //creates a map of the dead Ecal channels in the barrel and endcap
4238 //to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
4239 //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
4240 void
4241 OSUAnalysis::WriteDeadEcal (){
4242 double etaEcal, phiEcal;
4243 ifstream DeadEcalFile(deadEcalFile_);
4244 if(!DeadEcalFile) {
4245 clog << "Error: DeadEcalFile has not been found." << endl;
4246 return;
4247 }
4248 if(DeadEcalVec.size()!= 0){
4249 clog << "Error: DeadEcalVec has a nonzero size" << endl;
4250 return;
4251 }
4252 while(!DeadEcalFile.eof())
4253 {
4254 DeadEcalFile >> etaEcal >> phiEcal;
4255 DeadEcal newChan;
4256 newChan.etaEcal = etaEcal;
4257 newChan.phiEcal = phiEcal;
4258 DeadEcalVec.push_back(newChan);
4259 }
4260 if(DeadEcalVec.size() == 0) clog << "Warning: No dead Ecal channels have been found." << endl;
4261 }
4262
4263 //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
4264 int
4265 OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
4266 double deltaRLowest = 999;
4267 int value = 0;
4268 if (DeadEcalVec.size() == 0) WriteDeadEcal();
4269 for(vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
4270 double eta = ecal->etaEcal;
4271 double phi = ecal->phiEcal;
4272 double deltaRtemp = deltaR(eta, phi, track1->eta, track1->phi);
4273 if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
4274 }
4275 if (deltaRLowest<0.05) {value = 1;}
4276 else {value = 0;}
4277 return value;
4278 }
4279
4280 // Returns the smallest DeltaR between the object and any generated true particle in the event.
4281 template <class InputObject>
4282 double OSUAnalysis::getGenDeltaRLowest(InputObject object){
4283 double genDeltaRLowest = 999.;
4284 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
4285 double deltaRtemp = deltaR(mcparticle->eta, mcparticle->phi, object->eta, object->phi);
4286 if (deltaRtemp < genDeltaRLowest) genDeltaRLowest = deltaRtemp;
4287 }
4288 return genDeltaRLowest;
4289 }
4290
4291 double
4292 OSUAnalysis::applyFunction(string function, double value){
4293
4294 if(function == "abs") value = fabs(value);
4295 else if(function == "fabs") value = fabs(value);
4296 else if(function == "log10") value = log10(value);
4297 else if(function == "log") value = log10(value);
4298
4299 else if(function == "") value = value;
4300 else{clog << "WARNING: invalid function '" << function << "'\n";}
4301
4302 return value;
4303
4304 }
4305
4306
4307 template <class InputCollection>
4308 void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
4309
4310 if (verbose_>2) clog << " Beginning setObjectFlags for cut " << currentCutIndex << ": " << currentCut.name
4311 << ", inputType=" << inputType
4312 << endl;
4313 if (currentCut.inputCollection.find("pair")!=string::npos) {
4314 string obj1, obj2;
4315 getTwoObjs(currentCut.inputCollection, obj1, obj2);
4316 if (inputType==obj1 ||
4317 inputType==obj2) {
4318 // Do not add a cut to individualFlags or cumulativeFlags, if the cut is on a paired collection,
4319 // and the inputType is a member of the pair.
4320 // The cut will instead be applied when the setObjectFlags() is called for the paired collection.
4321 // For example, if currentCut.inputCollection==electron-muon pairs,
4322 // then the flags should not be set here when inputType==muons or inputType==electrons.
4323 return;
4324 }
4325 }
4326
4327 for (uint object = 0; object != inputCollection->size(); object++){
4328
4329 bool cutDecision = true;//object passes if this cut doesn't cut on that type of object
4330 bool plotDecision = true;
4331
4332 if(currentCut.inputCollection == inputType){
4333
4334 vector<bool> subcutDecisions;
4335 for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
4336 string stringValue = "";
4337 double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex), stringValue);
4338 if (stringValue == "") subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4339 else subcutDecisions.push_back(evaluateComparison(stringValue,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutStringValues.at(subcutIndex)));
4340
4341 }
4342 if(currentCut.numSubcuts == 1) cutDecision = subcutDecisions.at(0);
4343 else{
4344 bool tempDecision = true;
4345 for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
4346 if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
4347 tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
4348 else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
4349 tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
4350 }
4351 cutDecision = tempDecision;
4352 }
4353 //invert the cut if this cut is a veto
4354 if(currentCut.isVeto) cutDecision = !cutDecision;
4355 plotDecision = cutDecision;
4356 }
4357
4358 individualFlags.at(inputType).at(currentCutIndex).push_back(make_pair(cutDecision,plotDecision));
4359
4360 //set flags for objects that pass this cut AND all the previous cuts
4361 bool previousCumulativeCutFlag = true;
4362 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4363 if(previousCumulativeCutFlag && individualFlags.at(inputType).at(previousCutIndex).at(object).first) previousCumulativeCutFlag = true;
4364 else{ previousCumulativeCutFlag = false; break;}
4365 }
4366 previousCumulativeCutFlag = previousCumulativeCutFlag && cutDecision;
4367 bool previousCumulativePlotFlag = true;
4368 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4369 if(previousCumulativePlotFlag && individualFlags.at(inputType).at(previousCutIndex).at(object).second) previousCumulativePlotFlag = true;
4370 else{ previousCumulativePlotFlag = false; break;}
4371 }
4372 previousCumulativePlotFlag = previousCumulativePlotFlag && plotDecision;
4373
4374 cumulativeFlags.at(inputType).at(currentCutIndex).push_back(make_pair(previousCumulativeCutFlag,previousCumulativePlotFlag));
4375
4376 }
4377
4378 } // end void OSUAnalysis::setObjectFlags
4379
4380
4381 template <class InputCollection1, class InputCollection2>
4382 void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags,
4383 InputCollection1 inputCollection1, InputCollection2 inputCollection2, string inputType){
4384 // This function sets the flags for the paired object collection.
4385 // If the cut is applying on the given paired object collection, then the flags for the single object collections are also set.
4386 // If not, then the flags for the paired object collection are taken as the AND of the flags for each single object collection.
4387
4388 if (verbose_>2) clog << " Beginning setObjectFlags for cut=" << currentCut.name
4389 << ", inputType=" << inputType
4390 << endl;
4391
4392 bool sameObjects = false;
4393 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true; // FIXME: is sameObjects just the not of isTwoTypesOfObject? If so, it's redundant.
4394
4395 // Get the strings for the two objects that make up the pair.
4396 string obj1Type, obj2Type;
4397 getTwoObjs(inputType, obj1Type, obj2Type);
4398 bool isTwoTypesOfObject = true;
4399 if (obj1Type==obj2Type) isTwoTypesOfObject = false;
4400
4401 // Initialize the flags for individual objects to all be false, if the cut is on the pair.
4402 // Set them to true later, if any paired object passes (in which case both of its constituents should pass).
4403 if (currentCut.inputCollection == inputType) {
4404 for (uint object1 = 0; object1 != inputCollection1->size(); object1++) {
4405 individualFlags.at(obj1Type).at(currentCutIndex).push_back(make_pair(false,false));
4406 cumulativeFlags.at(obj1Type).at(currentCutIndex).push_back(make_pair(false,false));
4407 }
4408 if (isTwoTypesOfObject) { // Only initialize the second object if it is different from the first.
4409 for (uint object2 = 0; object2 != inputCollection2->size(); object2++) {
4410 individualFlags.at(obj2Type).at(currentCutIndex).push_back(make_pair(false,false));
4411 cumulativeFlags.at(obj2Type).at(currentCutIndex).push_back(make_pair(false,false));
4412 }
4413 }
4414 }
4415
4416 int counter = 0;
4417
4418 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
4419 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
4420
4421 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
4422
4423
4424 bool cutDecision = true;//object passes if this cut doesn't cut on that type of object
4425 bool plotDecision = true;
4426
4427 // Determine whether each pair passes the cut, only if inputCollection is the same as the inputType.
4428 if(currentCut.inputCollection == inputType){
4429
4430 vector<bool> subcutDecisions;
4431 for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
4432 string stringValue = "";
4433 double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex), stringValue);
4434 if (verbose_>1) clog << currentCut.variables.at(subcutIndex) << " = " << value
4435 << endl;
4436 if (stringValue == "") subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4437 else subcutDecisions.push_back(evaluateComparison(stringValue,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutStringValues.at(subcutIndex)));
4438 }
4439
4440 if(currentCut.numSubcuts == 1) cutDecision = subcutDecisions.at(0);
4441 else{
4442 bool tempDecision = subcutDecisions.at(0);
4443 for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
4444 if(currentCut.logicalOperators.at(subcutIndex-1) == "&" || currentCut.logicalOperators.at(subcutIndex-1) == "&&")
4445 tempDecision = tempDecision && subcutDecisions.at(subcutIndex);
4446 else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
4447 tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
4448 }
4449 cutDecision = tempDecision;
4450 }
4451 //invert the cut if this cut is a veto
4452 if (currentCut.isVeto) cutDecision = !cutDecision;
4453 plotDecision = cutDecision;
4454
4455 if (verbose_>1) clog << " cutDecision = " << cutDecision
4456 << "; for currentCut.inputCollection = " << currentCut.inputCollection
4457 << "; object1 (" << obj1Type << ") = " << object1
4458 << "; object2 (" << obj2Type << ") = " << object2
4459 << endl;
4460
4461 if (cutDecision) { // only set the flags for the individual objects if the pair object is being cut on
4462 individualFlags.at(obj1Type).at(currentCutIndex).at(object1).first = true;
4463 individualFlags.at(obj2Type).at(currentCutIndex).at(object2).first = true;
4464 }
4465 if (plotDecision) { // only set the flags for the individual objects if the pair object is being cut on
4466 individualFlags.at(obj1Type).at(currentCutIndex).at(object1).second = true;
4467 individualFlags.at(obj2Type).at(currentCutIndex).at(object2).second = true;
4468 }
4469
4470
4471 } // if(currentCut.inputCollection == inputType){
4472
4473 // The individualFlags will be true if the inputCollection is not the same as the inputType.
4474 // They are also independent of the previous flags on the single objects.
4475 individualFlags.at(inputType).at(currentCutIndex).push_back(make_pair(cutDecision,plotDecision));
4476
4477
4478
4479 // ************************************
4480 // Determine cumulative flags
4481 // ************************************
4482 // determine whether this paired object passes this cut AND all previous cuts
4483 bool previousCumulativeCutFlag = true;
4484 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4485 if(previousCumulativeCutFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter).first) previousCumulativeCutFlag = true;
4486 else{ previousCumulativeCutFlag = false; break;}
4487 }
4488 previousCumulativeCutFlag = previousCumulativeCutFlag && cutDecision;
4489
4490 bool previousCumulativePlotFlag = true;
4491 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4492 if(previousCumulativePlotFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter).second) previousCumulativePlotFlag = true;
4493 else{ previousCumulativePlotFlag = false; break;}
4494 }
4495 previousCumulativePlotFlag = previousCumulativePlotFlag && plotDecision;
4496
4497 // Get the index for the flags of each of the single objects in the pair. Usually this is the index of the previous cut, i.e., currentCutIndex-1.
4498 int cutIdxFlagsObj1 = max(int(currentCutIndex-1),0);
4499 int cutIdxFlagsObj2 = max(int(currentCutIndex-1),0);
4500 // If the inputCollection of the cut is not equal to the inputType but the inputCollection includes objects that are contained in inputType, then use the currentCutIndex for those collections.
4501 // For example, if the inputType is jet-jet pairs, and the inputCollection is track-jet pairs, then use currentCutIndex for cutIdxFlagsObj{1,2}, i.e., for both jets.
4502 // For example, if the inputType is jets, and the inputCollection is track-jet pairs, then use currentCutIndex for cutIdxFlagsObj2 (jets) and currentCutIndex-1 for cutIdxFlagsObj1 (tracks).
4503 if (currentCut.inputCollection != inputType) {
4504 if (currentCut.inputCollection.find(obj1Type)!=string::npos) cutIdxFlagsObj1 = currentCutIndex;
4505 if (currentCut.inputCollection.find(obj2Type)!=string::npos) cutIdxFlagsObj2 = currentCutIndex;
4506 }
4507 flagPair flags1 = cumulativeFlags.at(obj1Type).at(cutIdxFlagsObj1); // flag for input collection 1
4508 flagPair flags2 = cumulativeFlags.at(obj2Type).at(cutIdxFlagsObj2); // flag for input collection 2
4509
4510 // The cumulative flag is only true if the paired object cumulative flag is true, and if the single object cumulative flags are true.
4511 bool currentCumulativeCutFlag = true;
4512 bool currentCumulativePlotFlag = true;
4513 if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag;
4514 else if(flags1.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag && flags2.at(object2).first;
4515 else if(flags2.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag && flags1.at(object1).first;
4516 else currentCumulativeCutFlag = previousCumulativeCutFlag && flags1.at(object1).first && flags2.at(object2).first;
4517
4518 if(flags1.size() == 0 && flags2.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag;
4519 else if(flags1.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag && flags2.at(object2).second;
4520 else if(flags2.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag && flags1.at(object1).second;
4521 else currentCumulativePlotFlag = previousCumulativePlotFlag && flags1.at(object1).first && flags2.at(object2).second;
4522
4523 cumulativeFlags.at(inputType).at(currentCutIndex).push_back(make_pair(currentCumulativeCutFlag,currentCumulativePlotFlag)); // Set the flag for the paired object
4524
4525
4526 if (currentCumulativeCutFlag && currentCut.inputCollection == inputType) { // Set the flags for the individual objects if the paired object is being cut on.
4527 cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1).first = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1, "cut");
4528 cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2).first = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "cut");
4529 cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1).second = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1, "plot");
4530 cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2).second = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "plot");
4531
4532 if (verbose_>1) clog << " previousCumulativeCutFlag for object1 = " << getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object1, "cut") << endl;
4533 if (verbose_>1) clog << " previousCumulativeCutFlag for object2 = " << getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "cut") << endl;
4534
4535 }
4536
4537 counter++;
4538
4539 } // end for (uint object2 = 0; object2 != inputCollection2->size(); object2++)
4540 } // end for (uint object1 = 0; object1 != inputCollection1->size(); object1++)
4541
4542 } // end void OSUAnalysis::setObjectFlags
4543
4544
4545 bool OSUAnalysis::getPreviousCumulativeFlags(uint currentCutIndex, flagMap &individualFlags, string obj1Type, uint object1, string flagType) {
4546 // Return true iff for the collection obj1Type, the element with index object1 has individal flags set to true for
4547 // all cuts up to currentCutIndex
4548 bool previousCumulativeFlag = true;
4549 for (uint previousCutIndex = 0; previousCutIndex < currentCutIndex; previousCutIndex++) {
4550 bool tempFlag = false;
4551 if(flagType == "cut") tempFlag = individualFlags.at(obj1Type).at(previousCutIndex).at(object1).first;
4552 else if(flagType == "plot") tempFlag = individualFlags.at(obj1Type).at(previousCutIndex).at(object1).second;
4553
4554 if (previousCumulativeFlag && tempFlag) previousCumulativeFlag = true;
4555 else {
4556 previousCumulativeFlag = false; break;
4557 }
4558 }
4559 return previousCumulativeFlag;
4560 }
4561
4562
4563
4564
4565 template <class InputCollection>
4566 void OSUAnalysis::assignTreeBranch(BranchSpecs parameters, InputCollection inputCollection, flagPair flags){
4567 // This function is similar to fill1DHistogram(), but instead of filling a histogram it assigns a value to a variable for the BNTree
4568
4569 if (BNTreeBranchVals_.count(parameters.name)==0) clog << "Error[assignTreeBranch]: trying to assign value to " << parameters.name << " that does not have a branch set up. Will likely seg fault." << endl;
4570 for (uint object = 0; object != inputCollection->size(); object++) {
4571
4572 if (!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4573
4574 string inputVariable = parameters.inputVariable;
4575 string function = "";
4576 string stringValue = "";
4577 double value = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4578 BNTreeBranchVals_.at(parameters.name).push_back(value);
4579
4580 }
4581 }
4582
4583
4584 template <class InputCollection>
4585 void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection, flagPair flags, double scaleFactor){
4586
4587 if (verbose_>2) clog << " Filling histogram for " << parameters.name << endl;
4588
4589 for (uint object = 0; object != inputCollection->size(); object++){
4590
4591 if(!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4592
4593 string currentString = parameters.inputVariables.at(0);
4594 string inputVariable = "";
4595 string function = "";
4596 if(currentString.find("(")==string::npos){
4597 inputVariable = currentString;// variable to cut on
4598 }
4599 else{
4600 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4601 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4602 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4603 }
4604
4605 string stringValue = "";
4606 double value = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4607 histo->Fill(value,scaleFactor);
4608
4609 if (printEventInfo_) {
4610 // Write information about event to screen, for testing purposes.
4611 clog << " Info for event: value for histogram " << histo->GetName() << ": " << value << " (object number " << object << ")" << endl;
4612 }
4613
4614 }
4615 }
4616
4617 template <class InputCollection1, class InputCollection2>
4618 void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2,
4619 flagPair pairFlags, double scaleFactor){
4620
4621 bool sameObjects = false;
4622 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
4623
4624 int pairCounter = -1;
4625 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
4626 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
4627
4628 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
4629
4630 pairCounter++;
4631 if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter).second) continue;
4632
4633 string currentString = parameters.inputVariables.at(0);
4634 string inputVariable = "";
4635 string function = "";
4636 if(currentString.find("(")==string::npos){
4637 inputVariable = currentString;// variable to cut on
4638 }
4639 else{
4640 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4641 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4642 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4643 }
4644
4645 string stringValue = "";
4646 double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4647 histo->Fill(value,scaleFactor);
4648
4649 if (printEventInfo_) {
4650 // Write information about event to screen, for testing purposes.
4651 clog << " Info for event: value for histogram " << histo->GetName() << ": " << value
4652 << " (object1 number " << object1 << "), "
4653 << " (object2 number " << object2 << ")"
4654 << endl;
4655 }
4656
4657 }
4658 }
4659
4660 }
4661
4662
4663 template <class InputCollection>
4664 void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection, flagPair flags, double scaleFactor){
4665
4666 for (uint object = 0; object != inputCollection->size(); object++){
4667
4668 if(!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4669
4670 string stringValue = "";
4671 string currentString = parameters.inputVariables.at(0);
4672 string inputVariable = "";
4673 string function = "";
4674 if(currentString.find("(")==string::npos){
4675 inputVariable = currentString;// variable to cut on
4676 }
4677 else{
4678 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4679 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4680 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4681 }
4682 double valueX = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4683
4684 currentString = parameters.inputVariables.at(1);
4685 inputVariable = "";
4686 function = "";
4687 if(currentString.find("(")==string::npos){
4688 inputVariable = currentString;// variable to cut on
4689 }
4690 else{
4691 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4692 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4693 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4694 }
4695
4696 double valueY = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4697
4698 histo->Fill(valueX,valueY,scaleFactor);
4699
4700 }
4701
4702 }
4703
4704 template <class InputCollection1, class InputCollection2>
4705 void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2,
4706 flagPair pairFlags, double scaleFactor){
4707
4708 bool sameObjects = false;
4709 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
4710
4711 int pairCounter = -1;
4712 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
4713 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
4714
4715 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
4716
4717 pairCounter++;
4718
4719 if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter).second) continue;
4720
4721 string stringValue = "";
4722 string currentString = parameters.inputVariables.at(0);
4723 string inputVariable = "";
4724 string function = "";
4725 if(currentString.find("(")==string::npos){
4726 inputVariable = currentString;// variable to cut on
4727 }
4728 else{
4729 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4730 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4731 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4732 }
4733 double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4734
4735 currentString = parameters.inputVariables.at(1);
4736 inputVariable = "";
4737 function = "";
4738 if(currentString.find("(")==string::npos){
4739 inputVariable = currentString;// variable to cut on
4740 }
4741 else{
4742 function = currentString.substr(0,currentString.find("("));//function comes before the "("
4743 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4744 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4745 }
4746 double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4747
4748
4749 histo->Fill(valueX,valueY,scaleFactor);
4750
4751 }
4752 }
4753
4754 }
4755
4756
4757 template <class InputObject>
4758 int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
4759
4760 int bestMatchIndex = -1;
4761 double bestMatchDeltaR = 999;
4762
4763 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
4764
4765 double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
4766 if(currentDeltaR > 0.05) continue;
4767 // clog << setprecision(3) << setw(20)
4768 // << "\tcurrentParticle: eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
4769 // << setw(20)
4770 // << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
4771 // << setw(20)
4772 // << "\tdeltaR = " << currentDeltaR
4773 // << setprecision(1)
4774 // << setw(20)
4775 // << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
4776 // << setw(20)
4777 // << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
4778 // << setw(20)
4779 // << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
4780 if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
4781 bestMatchIndex = mcparticle - mcparticles->begin();
4782 bestMatchDeltaR = currentDeltaR;
4783 }
4784
4785 }
4786 // if(bestMatchDeltaR != 999) clog << "bestMatch: deltaR = " << bestMatchDeltaR << " id = " << mcparticles->at(bestMatchIndex).id << " motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
4787 // else clog << "no match found..." << endl;
4788 return bestMatchIndex;
4789
4790 }
4791
4792
4793 int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
4794
4795 int bestMatchIndex = -1;
4796 double bestMatchDeltaR = 999;
4797
4798 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
4799
4800 if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
4801
4802 double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
4803 if(currentDeltaR > 0.05) continue;
4804
4805 if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
4806 bestMatchIndex = mcparticle - mcparticles->begin();
4807 bestMatchDeltaR = currentDeltaR;
4808 }
4809
4810 }
4811
4812 return bestMatchIndex;
4813 }
4814
4815
4816 // bin particle type
4817 // --- -------------
4818 // 0 unmatched
4819 // 1 u
4820 // 2 d
4821 // 3 s
4822 // 4 c
4823 // 5 b
4824 // 6 t
4825 // 7 e
4826 // 8 mu
4827 // 9 tau
4828 // 10 nu
4829 // 11 g
4830 // 12 gamma
4831 // 13 Z
4832 // 14 W
4833 // 15 light meson
4834 // 16 K meson
4835 // 17 D meson
4836 // 18 B meson
4837 // 19 light baryon
4838 // 20 strange baryon
4839 // 21 charm baryon
4840 // 22 bottom baryon
4841 // 23 QCD string
4842 // 24 other
4843
4844 int OSUAnalysis::getPdgIdBinValue(int pdgId){
4845
4846 int binValue = -999;
4847 int absPdgId = fabs(pdgId);
4848 if(pdgId == -1) binValue = 0;
4849 else if(absPdgId <= 6 ) binValue = absPdgId;
4850 else if(absPdgId == 11 ) binValue = 7;
4851 else if(absPdgId == 13 ) binValue = 8;
4852 else if(absPdgId == 15 ) binValue = 9;
4853 else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
4854 else if(absPdgId == 21 ) binValue = 11;
4855 else if(absPdgId == 22 ) binValue = 12;
4856 else if(absPdgId == 23 ) binValue = 13;
4857 else if(absPdgId == 24 ) binValue = 14;
4858 else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130 ) binValue = 15;
4859 else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400) ) binValue = 16;
4860 else if(absPdgId > 400 && absPdgId < 500 ) binValue = 17;
4861 else if(absPdgId > 500 && absPdgId < 600 ) binValue = 18;
4862 else if(absPdgId > 1000 && absPdgId < 3000 ) binValue = 19;
4863 else if(absPdgId > 3000 && absPdgId < 4000 ) binValue = 20;
4864 else if(absPdgId > 4000 && absPdgId < 5000 ) binValue = 21;
4865 else if(absPdgId > 5000 && absPdgId < 6000 ) binValue = 22;
4866 else if(absPdgId == 92 ) binValue = 23;
4867
4868 else binValue = 24;
4869
4870 return binValue;
4871
4872 }
4873
4874 const BNprimaryvertex *
4875 OSUAnalysis::chosenVertex ()
4876 {
4877 const BNprimaryvertex *chosenVertex = 0;
4878 if(cumulativeFlags.find ("primaryvertexs") != cumulativeFlags.end ()){
4879 flagPair vertexFlags;
4880 for (int i = cumulativeFlags.at("primaryvertexs").size() - 1; i >= 0; i--){
4881 if (cumulativeFlags.at("primaryvertexs").at(i).size()){
4882 vertexFlags = cumulativeFlags.at("primaryvertexs").at(i);
4883 break;
4884 }
4885 }
4886 for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
4887 if(!vertexFlags.at(vertexIndex).first) continue;
4888 chosenVertex = & primaryvertexs->at(vertexIndex);
4889 break;
4890 }
4891 }
4892 else if (find (objectsToGet.begin (), objectsToGet.end (), "primaryvertexs") != objectsToGet.end ())
4893 chosenVertex = & primaryvertexs->at (0);
4894
4895 return chosenVertex;
4896 }
4897
4898 const BNmet *
4899 OSUAnalysis::chosenMET ()
4900 {
4901 const BNmet *chosenMET = 0;
4902 if(cumulativeFlags.find ("mets") != cumulativeFlags.end ()){
4903 flagPair metFlags;
4904 for (int i = cumulativeFlags.at("mets").size() - 1; i >= 0; i--){
4905 if (cumulativeFlags.at("mets").at(i).size()){
4906 metFlags = cumulativeFlags.at("mets").at(i);
4907 break;
4908 }
4909 }
4910 for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
4911 if(!metFlags.at(metIndex).first) continue;
4912 chosenMET = & mets->at(metIndex);
4913 break;
4914 }
4915 }
4916 else if (find (objectsToGet.begin (), objectsToGet.end (), "mets") != objectsToGet.end ())
4917 chosenMET = & mets->at (0);
4918
4919 return chosenMET;
4920 }
4921
4922 const BNelectron *
4923 OSUAnalysis::chosenElectron ()
4924 {
4925 const BNelectron *chosenElectron = 0;
4926 if(cumulativeFlags.find ("electrons") != cumulativeFlags.end ()){
4927 flagPair electronFlags;
4928 for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
4929 if (cumulativeFlags.at("electrons").at(i).size()){
4930 electronFlags = cumulativeFlags.at("electrons").at(i);
4931 break;
4932 }
4933 }
4934 for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
4935 if(!electronFlags.at(electronIndex).first) continue;
4936 chosenElectron = & electrons->at(electronIndex);
4937 break;
4938 }
4939 }
4940 else if (find (objectsToGet.begin (), objectsToGet.end (), "electrons") != objectsToGet.end ())
4941 chosenElectron = & electrons->at (0);
4942
4943 return chosenElectron;
4944 }
4945
4946
4947 const BNmuon *
4948 OSUAnalysis::chosenMuon ()
4949 {
4950 const BNmuon *chosenMuon = 0;
4951 if(cumulativeFlags.find ("muons") != cumulativeFlags.end ()){
4952 flagPair muonFlags;
4953 for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
4954 if (cumulativeFlags.at("muons").at(i).size()){
4955 muonFlags = cumulativeFlags.at("muons").at(i);
4956 break;
4957 }
4958 }
4959 for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
4960 if(!muonFlags.at(muonIndex).first) continue;
4961 chosenMuon = & muons->at(muonIndex);
4962 break;
4963 }
4964 }
4965 else if (find (objectsToGet.begin (), objectsToGet.end (), "muons") != objectsToGet.end ())
4966 chosenMuon = & muons->at (0);
4967
4968 return chosenMuon;
4969 }
4970
4971 double
4972 OSUAnalysis::chosenHT ()
4973 {
4974 double chosenHT = 0.0;
4975 if(cumulativeFlags.find ("jets") != cumulativeFlags.end ()){
4976 flagPair jetFlags;
4977 for (int i = cumulativeFlags.at("jets").size() - 1; i >= 0; i--){
4978 if (cumulativeFlags.at("jets").at(i).size()){
4979 jetFlags = cumulativeFlags.at("jets").at(i);
4980 break;
4981 }
4982 }
4983 for (uint jetIndex = 0; jetIndex != jetFlags.size(); jetIndex++){
4984 if(!jetFlags.at(jetIndex).first) continue;
4985 chosenHT += jets->at(jetIndex).pt;
4986 }
4987 }
4988
4989 return chosenHT;
4990 }
4991
4992 pair<const BNmuon *, const BNmuon *>
4993 OSUAnalysis::leadMuonPair ()
4994 {
4995 pair<const BNmuon *, const BNmuon *> leadMuonPair;
4996 leadMuonPair.first = leadMuonPair.second = 0;
4997
4998 if(cumulativeFlags.find ("muons") != cumulativeFlags.end ()){
4999 flagPair muonFlags;
5000 for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
5001 if (cumulativeFlags.at("muons").at(i).size()){
5002 muonFlags = cumulativeFlags.at("muons").at(i);
5003 break;
5004 }
5005 }
5006 for (uint muonIndex0 = 0; muonIndex0 != muonFlags.size(); muonIndex0++){
5007 if(!muonFlags.at(muonIndex0).first) continue;
5008 for (uint muonIndex1 = muonIndex0 + 1; muonIndex1 < muonFlags.size(); muonIndex1++){
5009 if(!muonFlags.at(muonIndex1).first) continue;
5010 const BNmuon *mu0 = & muons->at(muonIndex0);
5011 const BNmuon *mu1 = & muons->at(muonIndex1);
5012 if(leadMuonPair.first == 0 || leadMuonPair.second == 0){
5013 leadMuonPair.first = mu0;
5014 leadMuonPair.second = mu1;
5015 }
5016 else{
5017 TVector2 newPt0 (mu0->px, mu0->py),
5018 newPt1 (mu1->px, mu1->py),
5019 oldPt0 (leadMuonPair.first->px, leadMuonPair.first->py),
5020 oldPt1 (leadMuonPair.second->px, leadMuonPair.second->py);
5021 if(newPt0.Mod () + newPt1.Mod () > oldPt0.Mod() + oldPt1.Mod ())
5022 {
5023 leadMuonPair.first = mu0;
5024 leadMuonPair.second = mu1;
5025 }
5026 }
5027 }
5028 }
5029 }
5030
5031 return leadMuonPair;
5032 }
5033
5034 pair<const BNelectron *, const BNelectron *>
5035 OSUAnalysis::leadElectronPair ()
5036 {
5037 pair<const BNelectron *, const BNelectron *> leadElectronPair;
5038 leadElectronPair.first = leadElectronPair.second = 0;
5039 if(cumulativeFlags.find ("electrons") != cumulativeFlags.end ()){
5040 flagPair electronFlags;
5041 for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
5042 if (cumulativeFlags.at("electrons").at(i).size()){
5043 electronFlags = cumulativeFlags.at("electrons").at(i);
5044 break;
5045 }
5046 }
5047 for (uint electronIndex0 = 0; electronIndex0 != electronFlags.size(); electronIndex0++){
5048 if(!electronFlags.at(electronIndex0).first) continue;
5049 for (uint electronIndex1 = electronIndex0 + 1; electronIndex1 < electronFlags.size(); electronIndex1++){
5050 if(!electronFlags.at(electronIndex1).first) continue;
5051 const BNelectron *el0 = & electrons->at(electronIndex0);
5052 const BNelectron *el1 = & electrons->at(electronIndex1);
5053 if(leadElectronPair.first == 0 || leadElectronPair.second == 0){
5054 leadElectronPair.first = el0;
5055 leadElectronPair.second = el1;
5056 }
5057 else{
5058 TVector2 newPt0 (el0->px, el0->py),
5059 newPt1 (el1->px, el1->py),
5060 oldPt0 (leadElectronPair.first->px, leadElectronPair.first->py),
5061 oldPt1 (leadElectronPair.second->px, leadElectronPair.second->py);
5062 if(newPt0.Mod () + newPt1.Mod () > oldPt0.Mod() + oldPt1.Mod ())
5063 {
5064 leadElectronPair.first = el0;
5065 leadElectronPair.second = el1;
5066 }
5067 }
5068 }
5069 }
5070 }
5071
5072 return leadElectronPair;
5073 }
5074
5075 DEFINE_FWK_MODULE(OSUAnalysis);