ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc
Revision: 1.39
Committed: Wed Mar 27 16:52:25 2013 UTC (12 years, 1 month ago) by jbrinson
Content type: text/plain
Branch: MAIN
Changes since 1.38: +153 -12 lines
Log Message:
Added d0wrtPV and dZwrtPV to BNtrack; Added several valuelookup functions for pairs of objects

File Contents

# Content
1 #include "OSUT3Analysis/AnaTools/plugins/OSUAnalysis.h"
2
3 OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
4 // Retrieve parameters from the configuration file.
5 jets_ (cfg.getParameter<edm::InputTag> ("jets")),
6 muons_ (cfg.getParameter<edm::InputTag> ("muons")),
7 electrons_ (cfg.getParameter<edm::InputTag> ("electrons")),
8 events_ (cfg.getParameter<edm::InputTag> ("events")),
9 taus_ (cfg.getParameter<edm::InputTag> ("taus")),
10 mets_ (cfg.getParameter<edm::InputTag> ("mets")),
11 tracks_ (cfg.getParameter<edm::InputTag> ("tracks")),
12 genjets_ (cfg.getParameter<edm::InputTag> ("genjets")),
13 mcparticles_ (cfg.getParameter<edm::InputTag> ("mcparticles")),
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 puFile_ (cfg.getParameter<std::string> ("puFile")),
20 deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")),
21
22 muonSFFile_ (cfg.getParameter<std::string> ("muonSFFile")),
23 dataPU_ (cfg.getParameter<std::string> ("dataPU")),
24 electronSFID_ (cfg.getParameter<std::string> ("electronSFID")),
25 muonSF_ (cfg.getParameter<std::string> ("muonSF")),
26 dataset_ (cfg.getParameter<std::string> ("dataset")),
27 datasetType_ (cfg.getParameter<std::string> ("datasetType")),
28 channels_ (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
29 histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
30 plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
31 doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting")),
32 printEventInfo_ (cfg.getParameter<bool> ("printEventInfo"))
33 {
34
35 TH1::SetDefaultSumw2 ();
36
37 //create pile-up reweighting object, if necessary
38 if(datasetType_ != "data") {
39 if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
40 // muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
41 // electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
42 }
43
44
45 // Construct Cutflow Objects. These store the results of cut decisions and
46 // handle filling cut flow histograms.
47 masterCutFlow_ = new CutFlow (fs_);
48 std::vector<TFileDirectory> directories;
49
50 //always get vertex collection so we can assign the primary vertex in the event
51 objectsToGet.push_back("primaryvertexs");
52
53 //always make the plot of number of primary verticex (to check pile-up reweighting)
54 objectsToPlot.push_back("primaryvertexs");
55
56 //always get the MC particles to do GEN-matching
57 objectsToGet.push_back("mcparticles");
58
59 //always get the event collection to do pile-up reweighting
60 objectsToGet.push_back("events");
61
62 //parse the histogram definitions
63 for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
64
65 string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
66 if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
67 if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
68 objectsToGet.push_back(tempInputCollection);
69 objectsToPlot.push_back(tempInputCollection);
70 objectsToCut.push_back(tempInputCollection);
71 }
72 else{//pair of objects, need to add them both to the things to objectsToGet
73 int dashIndex = tempInputCollection.find("-");
74 int spaceIndex = tempInputCollection.find(" ");
75 int secondWordLength = spaceIndex - dashIndex;
76 objectsToGet.push_back(tempInputCollection);
77 objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
78 objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
79 objectsToPlot.push_back(tempInputCollection);
80 objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
81 objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
82 objectsToCut.push_back(tempInputCollection);
83 objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
84 objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
85
86 }
87
88 vector<edm::ParameterSet> histogramList_ (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
89
90 for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
91
92 histogram tempHistogram;
93 tempHistogram.inputCollection = tempInputCollection;
94 tempHistogram.name = histogramList_.at(currentHistogram).getParameter<string>("name");
95 tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
96 tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
97 tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
98
99 histograms.push_back(tempHistogram);
100
101 }
102 }
103 //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
104 sort( objectsToPlot.begin(), objectsToPlot.end() );
105 objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
106
107
108
109 //add histograms with the gen-matched id, mother id, and grandmother id
110 for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
111
112 string currentObject = objectsToPlot.at(currentObjectIndex);
113 if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
114
115 histogram tempIdHisto;
116 histogram tempMomIdHisto;
117 histogram tempGmaIdHisto;
118 histogram tempIdVsMomIdHisto;
119
120 tempIdHisto.inputCollection = currentObject;
121 tempMomIdHisto.inputCollection = currentObject;
122 tempGmaIdHisto.inputCollection = currentObject;
123 tempIdVsMomIdHisto.inputCollection = currentObject;
124
125 currentObject = currentObject.substr(0, currentObject.size()-1);
126 tempIdHisto.name = currentObject+"GenMatchId";
127 tempMomIdHisto.name = currentObject+"GenMatchMotherId";
128 tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
129 tempIdVsMomIdHisto.name = currentObject+"GenMatchIdVsMotherId";
130
131 currentObject.at(0) = toupper(currentObject.at(0));
132 tempIdHisto.title = currentObject+" Gen-matched Particle";
133 tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
134 tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
135 tempIdVsMomIdHisto.title = currentObject+" Gen-matched Particle's Mother vs. Particle;Particle;Mother";
136
137 int maxNum = 24;
138 vector<double> binVector;
139 binVector.push_back(maxNum);
140 binVector.push_back(0);
141 binVector.push_back(maxNum);
142
143 tempIdHisto.bins = binVector;
144 tempIdHisto.inputVariables.push_back("genMatchedId");
145 tempMomIdHisto.bins = binVector;
146 tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
147 tempGmaIdHisto.bins = binVector;
148 tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
149 binVector.push_back(maxNum);
150 binVector.push_back(0);
151 binVector.push_back(maxNum);
152 tempIdVsMomIdHisto.bins = binVector;
153 tempIdVsMomIdHisto.inputVariables.push_back("genMatchedId");
154 tempIdVsMomIdHisto.inputVariables.push_back("genMatchedMotherIdReverse");
155
156 histograms.push_back(tempIdHisto);
157 histograms.push_back(tempMomIdHisto);
158 histograms.push_back(tempGmaIdHisto);
159 histograms.push_back(tempIdVsMomIdHisto);
160 }
161
162
163 channel tempChannel;
164 //loop over all channels (event selections)
165 for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
166
167 //get name of channel
168 string channelName (channels_.at(currentChannel).getParameter<string>("name"));
169 tempChannel.name = channelName;
170 TString channelLabel = channelName;
171
172 //set triggers for this channel
173 vector<string> triggerNames;
174 triggerNames.clear();
175 tempChannel.triggers.clear();
176 if(channels_.at(currentChannel).exists("triggers")){
177 triggerNames = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
178 for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
179 tempChannel.triggers.push_back(triggerNames.at(trigger));
180 objectsToGet.push_back("triggers");
181 }
182
183
184
185
186 //create cutFlow for this channel
187 cutFlows_.push_back (new CutFlow (fs_, channelName));
188
189 //book a directory in the output file with the name of the channel
190 TFileDirectory subDir = fs_->mkdir( channelName );
191 directories.push_back(subDir);
192
193 std::map<std::string, TH1D*> oneDhistoMap;
194 oneDHists_.push_back(oneDhistoMap);
195 std::map<std::string, TH2D*> twoDhistoMap;
196 twoDHists_.push_back(twoDhistoMap);
197
198
199
200 //book all histograms included in the configuration
201 for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
202 histogram currentHistogram = histograms.at(currentHistogramIndex);
203 int numBinsElements = currentHistogram.bins.size();
204 int numInputVariables = currentHistogram.inputVariables.size();
205
206 if(numBinsElements != 3 && numBinsElements !=6) {
207 std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
208 exit(0);
209 }
210 else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
211 std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
212 exit(0);
213 }
214 else if(numBinsElements == 3){
215 oneDHists_.at(currentChannel)[currentHistogram.name] = directories.at(currentChannel).make<TH1D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, currentHistogram.bins.at(0), currentHistogram.bins.at(1), currentHistogram.bins.at(2));
216 }
217 else if(numBinsElements == 6){
218 twoDHists_.at(currentChannel)[currentHistogram.name] = directories.at(currentChannel).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));
219
220 }
221
222
223 if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
224
225 // bin particle type
226 // --- -------------
227 // 0 unmatched
228 // 1 u
229 // 2 d
230 // 3 s
231 // 4 c
232 // 5 b
233 // 6 t
234 // 7 e
235 // 8 mu
236 // 9 tau
237 // 10 nu
238 // 11 g
239 // 12 gamma
240 // 13 Z
241 // 14 W
242 // 15 light meson
243 // 16 K meson
244 // 17 D meson
245 // 18 B meson
246 // 19 light baryon
247 // 20 strange baryon
248 // 21 charm baryon
249 // 22 bottom baryon
250 // 23 other
251
252 vector<TString> labelArray;
253 labelArray.push_back("unmatched");
254 labelArray.push_back("u");
255 labelArray.push_back("d");
256 labelArray.push_back("s");
257 labelArray.push_back("c");
258 labelArray.push_back("b");
259 labelArray.push_back("t");
260 labelArray.push_back("e");
261 labelArray.push_back("#mu");
262 labelArray.push_back("#tau");
263 labelArray.push_back("#nu");
264 labelArray.push_back("g");
265 labelArray.push_back("#gamma");
266 labelArray.push_back("Z");
267 labelArray.push_back("W");
268 labelArray.push_back("light meson");
269 labelArray.push_back("K meson");
270 labelArray.push_back("D meson");
271 labelArray.push_back("B meson");
272 labelArray.push_back("light baryon");
273 labelArray.push_back("strange baryon");
274 labelArray.push_back("charm baryon");
275 labelArray.push_back("bottom baryon");
276 labelArray.push_back("other");
277
278 for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
279 if(currentHistogram.name.find("GenMatchIdVsMotherId")==std::string::npos) {
280 oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
281 }
282 else {
283 twoDHists_.at(currentChannel)[currentHistogram.name]->GetYaxis()->SetBinLabel(bin+1,labelArray.at(currentHistogram.bins.at(0)-bin-1));
284 twoDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
285 }
286 }
287 if(currentHistogram.name.find("GenMatchIdVsMotherId")!=std::string::npos) {
288 twoDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->CenterTitle();
289 twoDHists_.at(currentChannel)[currentHistogram.name]->GetYaxis()->CenterTitle();
290 }
291
292 }
293
294 //book a histogram for the number of each object type to be plotted
295
296 for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
297 string currentObject = objectsToPlot.at(currentObjectIndex);
298 int maxNum = 10;
299 if(currentObject == "mcparticles") maxNum = 50;
300 else if(currentObject == "primaryvertexs") maxNum = 50;
301 else if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
302 else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
303 else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
304
305 else if(currentObject == "electron-track pairs") currentObject = "electronTrackPairs";
306 else if(currentObject == "muon-track pairs") currentObject = "muonTrackPairs";
307 else if(currentObject == "muon-tau pairs") currentObject = "muonTauPairs";
308 else if(currentObject == "tau-tau pairs") currentObject = "ditauPairs";
309
310 currentObject.at(0) = toupper(currentObject.at(0));
311 string histoName = "num" + currentObject;
312
313 if(histoName == "numPrimaryvertexs"){
314 string newHistoName = histoName + "BeforePileupCorrection";
315 oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
316 newHistoName = histoName + "AfterPileupCorrection";
317 oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
318 }
319 else
320 oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
321 }
322
323
324
325
326 //get list of cuts for this channel
327 vector<edm::ParameterSet> cuts_ (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
328
329 //loop over and parse all cuts
330 for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
331 cut tempCut;
332 //store input collection for cut
333 string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
334 tempCut.inputCollection = tempInputCollection;
335 if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
336 objectsToGet.push_back(tempInputCollection);
337 objectsToCut.push_back(tempInputCollection);
338 }
339 else{//pair of objects, need to add them both to the things to objectsToGet
340 int dashIndex = tempInputCollection.find("-");
341 int spaceIndex = tempInputCollection.find(" ");
342 int secondWordLength = spaceIndex - dashIndex;
343 objectsToGet.push_back(tempInputCollection);
344 objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
345 objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
346 objectsToCut.push_back(tempInputCollection);
347 objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
348 objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
349
350 }
351
352
353
354 //split cut string into parts and store them
355 string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
356 std::vector<string> cutStringVector = splitString(cutString);
357 if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
358 std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
359 exit(0);
360 }
361 tempCut.numSubcuts = (cutStringVector.size()+1)/4;
362 for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
363 int indexOffset = 4 * subcutIndex;
364 string currentVariableString = cutStringVector.at(indexOffset);
365 if(currentVariableString.find("(")==std::string::npos){
366 tempCut.functions.push_back("");//no function was specified
367 tempCut.variables.push_back(currentVariableString);// variable to cut on
368 }
369 else{
370 tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
371 string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
372 tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
373 }
374 tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
375 tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
376 if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
377 }
378
379 //get number of objects required to pass cut for event to pass
380 string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
381 std::vector<string> numberRequiredVector = splitString(numberRequiredString);
382 if(numberRequiredVector.size()!=2){
383 std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
384 exit(0);
385 }
386
387 int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
388 tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
389 tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
390
391
392 string tempCutName;
393 if(cuts_.at(currentCut).exists("alias")){
394 tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
395 }
396 else{
397 //construct string for cutflow table
398 bool plural = numberRequiredInt != 1;
399 string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
400 string cutName = numberRequiredString + " " + collectionString + " with " + cutString;
401 tempCutName = cutName;
402 }
403 tempCut.name = tempCutName;
404
405
406 tempChannel.cuts.push_back(tempCut);
407
408
409
410 }//end loop over cuts
411
412 channels.push_back(tempChannel);
413 tempChannel.cuts.clear();
414
415 }//end loop over channels
416
417
418 //make unique vector of objects we need to get from the event
419 sort( objectsToGet.begin(), objectsToGet.end() );
420 objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
421 //make unique vector of objects we need to cut on
422 sort( objectsToCut.begin(), objectsToCut.end() );
423 objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
424
425
426 }
427
428 OSUAnalysis::~OSUAnalysis ()
429 {
430 // Destroying the CutFlow objects causes the cut flow numbers and time
431 // information to be printed to standard output.
432 for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
433 delete cutFlows_.at(currentChannel);
434 }
435 }
436
437 void
438 OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
439 {
440
441 // Retrieve necessary collections from the event.
442
443 if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
444 event.getByLabel (triggers_, triggers);
445
446 if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
447 event.getByLabel (jets_, jets);
448
449 if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
450 event.getByLabel (muons_, muons);
451
452 if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
453 event.getByLabel (electrons_, electrons);
454
455 if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
456 event.getByLabel (events_, events);
457
458 if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
459 event.getByLabel (taus_, taus);
460
461 if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
462 event.getByLabel (mets_, mets);
463
464 if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
465 event.getByLabel (tracks_, tracks);
466
467 if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
468 event.getByLabel (genjets_, genjets);
469
470 if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
471 event.getByLabel (mcparticles_, mcparticles);
472
473 if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
474 event.getByLabel (primaryvertexs_, primaryvertexs);
475
476 if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
477 event.getByLabel (bxlumis_, bxlumis);
478
479 if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
480 event.getByLabel (photons_, photons);
481
482 if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
483 event.getByLabel (superclusters_, superclusters);
484
485 //get pile-up event weight
486 double scaleFactor = 1.00;
487 if(doPileupReweighting_ && datasetType_ != "data")
488 scaleFactor = puWeight_->at (events->at (0).numTruePV);
489
490
491 //loop over all channels
492
493 for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
494 channel currentChannel = channels.at(currentChannelIndex);
495
496 flagMap individualFlags;
497 counterMap passingCounter;
498 cumulativeFlags.clear ();
499
500 bool triggerDecision = true;
501 if(currentChannel.triggers.size() != 0){ //triggers specified
502 triggerDecision = evaluateTriggers(currentChannel.triggers,triggers.product());
503 cutFlows_.at(currentChannelIndex)->at ("trigger") = triggerDecision;
504 }
505
506 //loop over all cuts
507
508
509
510 for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
511 cut currentCut = currentChannel.cuts.at(currentCutIndex);
512
513 for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
514 string currentObject = objectsToCut.at(currentObjectIndex);
515
516 //initialize maps to get ready to set cuts
517 individualFlags[currentObject].push_back (vector<bool> ());
518 cumulativeFlags[currentObject].push_back (vector<bool> ());
519
520 }
521 for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
522 string currentObject = objectsToCut.at(currentObjectIndex);
523
524 int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
525
526
527 if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
528 else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
529 else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
530 else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
531 else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
532 else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
533 else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
534 else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
535 else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
536 else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
537 else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
538 else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
539 else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
540
541
542 else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
543 cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
544 cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
545 "muon-muon pairs");
546 else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
547 cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
548 cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
549 "electron-electron pairs");
550 else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
551 cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
552 cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
553 "electron-muon pairs");
554 else if(currentObject == "electron-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),tracks.product(), \
555 cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
556 cumulativeFlags.at("tracks").at(flagsForPairCutsIndex), \
557 "electron-track pairs");
558 else if(currentObject == "muon-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),tracks.product(), \
559 cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
560 cumulativeFlags.at("tracks").at(flagsForPairCutsIndex), \
561 "muon-track pairs");
562 else if(currentObject == "muon-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),taus.product(), \
563 cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
564 cumulativeFlags.at("taus").at(flagsForPairCutsIndex), \
565 "muon-tau pairs");
566
567 else if(currentObject == "tau-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus .product(),taus.product(), \
568 cumulativeFlags.at("taus").at(flagsForPairCutsIndex), \
569 cumulativeFlags.at("taus").at(flagsForPairCutsIndex), \
570 "tau-tau pairs");
571
572
573 }
574
575
576
577 }//end loop over all cuts
578
579
580 //use cumulative flags to apply cuts at event level
581
582 bool eventPassedAllCuts = true;
583
584 //apply trigger (true if none were specified)
585 eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
586
587
588
589 for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
590
591 //loop over all objects and count how many passed the cumulative selection up to this point
592 cut currentCut = currentChannel.cuts.at(currentCutIndex);
593 int numberPassing = 0;
594
595 for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
596 if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
597 }
598
599 bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
600 cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
601
602 eventPassedAllCuts = eventPassedAllCuts && cutDecision;
603
604 }
605
606 cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
607
608
609
610 if(!eventPassedAllCuts)continue;
611
612 if (printEventInfo_) {
613 // Write information about event to screen, for testing purposes.
614 cout << "Event passed all cuts in channel " << currentChannel.name
615 << ": run=" << events->at(0).run
616 << " lumi=" << events->at(0).lumi
617 << " event=" << events->at(0).evt
618 << endl;
619 }
620
621
622 // if(datasetType_ != "data") {
623 // scaleFactor *= muonSFWeight_->at (chosenMuon ()->eta);
624 // scaleFactor *= electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
625 // }
626
627 //filling histograms
628 for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
629 histogram currentHistogram = histograms.at(histogramIndex);
630
631 if(currentHistogram.inputVariables.size() == 1){
632 TH1D* histo;
633 histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
634
635
636
637 if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
638 else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
639 else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
640 cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
641 cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
642 else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
643 else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
644 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
645 cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
646 else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
647 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
648 cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
649 else if(currentHistogram.inputCollection == "electron-track pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),tracks.product(), \
650 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("tracks").back(),
651 cumulativeFlags.at("electron-track pairs").back(),scaleFactor);
652 else if(currentHistogram.inputCollection == "muon-track pairs") fill1DHistogram(histo,currentHistogram, muons.product(),tracks.product(), \
653 cumulativeFlags.at("muons").back(),cumulativeFlags.at("tracks").back(),
654 cumulativeFlags.at("muon-track pairs").back(),scaleFactor);
655 else if(currentHistogram.inputCollection == "muon-tau pairs") fill1DHistogram(histo,currentHistogram, muons.product(),taus.product(), \
656 cumulativeFlags.at("muons").back(),cumulativeFlags.at("taus").back(),
657 cumulativeFlags.at("muon-tau pairs").back(),scaleFactor);
658
659 else if(currentHistogram.inputCollection == "tau-tau pairs") fill1DHistogram(histo,currentHistogram, taus.product(),taus.product(), \
660 cumulativeFlags.at("taus").back(),cumulativeFlags.at("taus").back(),
661 cumulativeFlags.at("tau-tau pairs").back(),scaleFactor);
662 else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
663 else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
664 else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
665 else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
666 else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
667 else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
668 else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
669 else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
670 else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
671 else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
672 }
673 else if(currentHistogram.inputVariables.size() == 2){
674 TH2D* histo;
675 histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
676
677
678
679 if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
680 else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
681 else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
682 cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
683 cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
684 else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
685 else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
686 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
687 cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
688 else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
689 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
690 cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
691 else if(currentHistogram.inputCollection == "electron-track pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),tracks.product(), \
692 cumulativeFlags.at("electrons").back(),cumulativeFlags.at("tracks").back(), \
693 cumulativeFlags.at("electron-track pairs").back(),scaleFactor);
694 else if(currentHistogram.inputCollection == "muon-track pairs") fill2DHistogram(histo,currentHistogram,muons.product(),tracks.product(), \
695 cumulativeFlags.at("muons").back(),cumulativeFlags.at("tracks").back(), \
696 cumulativeFlags.at("muon-track pairs").back(),scaleFactor);
697 else if(currentHistogram.inputCollection == "muon-tau pairs") fill2DHistogram(histo,currentHistogram,muons.product(),taus.product(), \
698 cumulativeFlags.at("muons").back(),cumulativeFlags.at("taus").back(), \
699 cumulativeFlags.at("muon-tau pairs").back(),scaleFactor);
700 else if(currentHistogram.inputCollection == "tau-tau pairs") fill2DHistogram(histo,currentHistogram,taus.product(),taus.product(), \
701 cumulativeFlags.at("taus").back(),cumulativeFlags.at("taus").back(), \
702 cumulativeFlags.at("tau-tau pairs").back(),scaleFactor);
703 else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
704 else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
705 else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
706 else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
707 else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
708 else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
709 else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
710 else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
711 else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
712 else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
713 }
714 }
715
716 //fills histograms with the sizes of collections
717 for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
718 string currentObject = objectsToPlot.at(currentObjectIndex);
719
720 string objectToPlot = "";
721
722 if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
723 else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
724 else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
725 else if(currentObject == "electron-track pairs") objectToPlot = "electronTrackPairs";
726 else if(currentObject == "muon-track pairs") objectToPlot = "muonTrackPairs";
727 else if(currentObject == "muon-tau pairs") objectToPlot = "muonTauPairs";
728 else if(currentObject == "tau-tau pairs") objectToPlot = "ditauPairs";
729 else objectToPlot = currentObject;
730 string tempCurrentObject = objectToPlot;
731 tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
732 string histoName = "num" + tempCurrentObject;
733
734
735
736
737 //set position of primary vertex in event, in order to calculate quantities relative to it
738 if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
739 vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
740 int numToPlot = 0;
741 for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
742 if(lastCutFlags.at(currentFlag)) numToPlot++;
743 }
744 if(objectToPlot == "primaryvertexs"){
745 oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
746 oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
747 }
748 else
749 oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,scaleFactor);
750 }
751 else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),scaleFactor);
752 else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),scaleFactor);
753 else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,scaleFactor);
754 else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),scaleFactor);
755 else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,scaleFactor);
756 else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),scaleFactor);
757 else if(objectToPlot == "electronTrackPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*tracks->size(),scaleFactor);
758 else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),scaleFactor);
759 else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),scaleFactor);
760 else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),scaleFactor);
761 else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),scaleFactor);
762 else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),scaleFactor);
763 else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),scaleFactor);
764 else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),scaleFactor);
765 else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),scaleFactor);
766 else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),scaleFactor);
767 else if(objectToPlot == "primaryvertexs"){
768 oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
769 oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
770 }
771
772 }
773
774
775
776
777 } //end loop over channel
778
779 masterCutFlow_->fillCutFlow(scaleFactor);
780
781
782
783 }
784
785
786 bool
787 OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
788
789
790 if(comparison == ">") return testValue > cutValue;
791 else if(comparison == ">=") return testValue >= cutValue;
792 else if(comparison == "<") return testValue < cutValue;
793 else if(comparison == "<=") return testValue <= cutValue;
794 else if(comparison == "==") return testValue == cutValue;
795 else if(comparison == "=") return testValue == cutValue;
796 else if(comparison == "!=") return testValue != cutValue;
797 else {std::cout << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
798
799 }
800
801 bool
802 OSUAnalysis::evaluateTriggers (vector<string> triggersToTest, const BNtriggerCollection* triggerCollection){
803
804 bool triggerDecision = false;
805
806 for (BNtriggerCollection::const_iterator trigger = triggerCollection->begin (); trigger != triggerCollection->end (); trigger++)
807 {
808 if(trigger->pass != 1) continue;
809 for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++)
810 {
811 if(trigger->name.find(triggersToTest.at(triggerName))!=std::string::npos){
812 triggerDecision = true;
813 }
814 }
815 }
816 return triggerDecision;
817
818 }
819
820 std::vector<std::string>
821 OSUAnalysis::splitString (string inputString){
822
823 std::stringstream stringStream(inputString);
824 std::istream_iterator<std::string> begin(stringStream);
825 std::istream_iterator<std::string> end;
826 std::vector<std::string> stringVector(begin, end);
827 return stringVector;
828
829 }
830
831 double
832 OSUAnalysis::valueLookup (const BNjet* object, string variable, string function){
833
834 double value = 0.0;
835 if(variable == "energy") value = object->energy;
836 else if(variable == "et") value = object->et;
837 else if(variable == "pt") value = object->pt;
838 else if(variable == "px") value = object->px;
839 else if(variable == "py") value = object->py;
840 else if(variable == "pz") value = object->pz;
841 else if(variable == "phi") value = object->phi;
842 else if(variable == "eta") value = object->eta;
843 else if(variable == "theta") value = object->theta;
844 else if(variable == "Upt") value = object->Upt;
845 else if(variable == "Uenergy") value = object->Uenergy;
846 else if(variable == "L2pt") value = object->L2pt;
847 else if(variable == "L2L3pt") value = object->L2L3pt;
848 else if(variable == "L2L3respt") value = object->L2L3respt;
849 else if(variable == "respt") value = object->respt;
850 else if(variable == "EMfrac") value = object->EMfrac;
851 else if(variable == "Hadfrac") value = object->Hadfrac;
852 else if(variable == "charge") value = object->charge;
853 else if(variable == "mass") value = object->mass;
854 else if(variable == "area") value = object->area;
855 else if(variable == "fHPD") value = object->fHPD;
856 else if(variable == "approximatefHPD") value = object->approximatefHPD;
857 else if(variable == "genPartonET") value = object->genPartonET;
858 else if(variable == "genPartonPT") value = object->genPartonPT;
859 else if(variable == "genPartonEta") value = object->genPartonEta;
860 else if(variable == "genPartonPhi") value = object->genPartonPhi;
861 else if(variable == "genJetET") value = object->genJetET;
862 else if(variable == "genJetPT") value = object->genJetPT;
863 else if(variable == "genJetEta") value = object->genJetEta;
864 else if(variable == "genJetPhi") value = object->genJetPhi;
865 else if(variable == "btagTChighPur") value = object->btagTChighPur;
866 else if(variable == "btagTChighEff") value = object->btagTChighEff;
867 else if(variable == "btagJetProb") value = object->btagJetProb;
868 else if(variable == "btagJetBProb") value = object->btagJetBProb;
869 else if(variable == "btagSoftEle") value = object->btagSoftEle;
870 else if(variable == "btagSoftMuon") value = object->btagSoftMuon;
871 else if(variable == "btagSoftMuonNoIP") value = object->btagSoftMuonNoIP;
872 else if(variable == "btagSecVertex") value = object->btagSecVertex;
873 else if(variable == "btagSecVertexHighEff") value = object->btagSecVertexHighEff;
874 else if(variable == "btagSecVertexHighPur") value = object->btagSecVertexHighPur;
875 else if(variable == "btagCombinedSecVertex") value = object->btagCombinedSecVertex;
876 else if(variable == "btagCombinedSecVertexMVA") value = object->btagCombinedSecVertexMVA;
877 else if(variable == "btagSoftMuonByPt") value = object->btagSoftMuonByPt;
878 else if(variable == "btagSoftMuonByIP3") value = object->btagSoftMuonByIP3;
879 else if(variable == "btagSoftElectronByPt") value = object->btagSoftElectronByPt;
880 else if(variable == "btagSoftElectronByIP3") value = object->btagSoftElectronByIP3;
881 else if(variable == "n90Hits") value = object->n90Hits;
882 else if(variable == "hitsInN90") value = object->hitsInN90;
883 else if(variable == "chargedHadronEnergyFraction") value = object->chargedHadronEnergyFraction;
884 else if(variable == "neutralHadronEnergyFraction") value = object->neutralHadronEnergyFraction;
885 else if(variable == "chargedEmEnergyFraction") value = object->chargedEmEnergyFraction;
886 else if(variable == "neutralEmEnergyFraction") value = object->neutralEmEnergyFraction;
887 else if(variable == "fLong") value = object->fLong;
888 else if(variable == "fShort") value = object->fShort;
889 else if(variable == "etaetaMoment") value = object->etaetaMoment;
890 else if(variable == "phiphiMoment") value = object->phiphiMoment;
891 else if(variable == "JESunc") value = object->JESunc;
892 else if(variable == "JECuncUp") value = object->JECuncUp;
893 else if(variable == "JECuncDown") value = object->JECuncDown;
894 else if(variable == "puJetMVA_full") value = object->puJetMVA_full;
895 else if(variable == "puJetMVA_simple") value = object->puJetMVA_simple;
896 else if(variable == "puJetMVA_cutbased") value = object->puJetMVA_cutbased;
897 else if(variable == "dZ") value = object->dZ;
898 else if(variable == "dR2Mean") value = object->dR2Mean;
899 else if(variable == "dRMean") value = object->dRMean;
900 else if(variable == "frac01") value = object->frac01;
901 else if(variable == "frac02") value = object->frac02;
902 else if(variable == "frac03") value = object->frac03;
903 else if(variable == "frac04") value = object->frac04;
904 else if(variable == "frac05") value = object->frac05;
905 else if(variable == "frac06") value = object->frac06;
906 else if(variable == "frac07") value = object->frac07;
907 else if(variable == "beta") value = object->beta;
908 else if(variable == "betaStar") value = object->betaStar;
909 else if(variable == "betaClassic") value = object->betaClassic;
910 else if(variable == "betaStarClassic") value = object->betaStarClassic;
911 else if(variable == "ptD") value = object->ptD;
912 else if(variable == "nvtx") value = object->nvtx;
913 else if(variable == "d0") value = object->d0;
914 else if(variable == "leadCandPt") value = object->leadCandPt;
915 else if(variable == "leadCandVx") value = object->leadCandVx;
916 else if(variable == "leadCandVy") value = object->leadCandVy;
917 else if(variable == "leadCandVz") value = object->leadCandVz;
918 else if(variable == "leadCandDistFromPV") value = object->leadCandDistFromPV;
919 else if(variable == "flavour") value = object->flavour;
920 else if(variable == "Nconst") value = object->Nconst;
921 else if(variable == "jetIDMinimal") value = object->jetIDMinimal;
922 else if(variable == "jetIDLooseAOD") value = object->jetIDLooseAOD;
923 else if(variable == "jetIDLoose") value = object->jetIDLoose;
924 else if(variable == "jetIDTight") value = object->jetIDTight;
925 else if(variable == "genPartonId") value = object->genPartonId;
926 else if(variable == "genPartonMotherId") value = object->genPartonMotherId;
927 else if(variable == "genPartonMother0Id") value = object->genPartonMother0Id;
928 else if(variable == "genPartonMother1Id") value = object->genPartonMother1Id;
929 else if(variable == "genPartonGrandMotherId") value = object->genPartonGrandMotherId;
930 else if(variable == "genPartonGrandMother00Id") value = object->genPartonGrandMother00Id;
931 else if(variable == "genPartonGrandMother01Id") value = object->genPartonGrandMother01Id;
932 else if(variable == "genPartonGrandMother10Id") value = object->genPartonGrandMother10Id;
933 else if(variable == "genPartonGrandMother11Id") value = object->genPartonGrandMother11Id;
934 else if(variable == "chargedMultiplicity") value = object->chargedMultiplicity;
935 else if(variable == "neutralMultiplicity") value = object->neutralMultiplicity;
936 else if(variable == "nconstituents") value = object->nconstituents;
937 else if(variable == "nHit") value = object->nHit;
938 else if(variable == "puJetId_full") value = object->puJetId_full;
939 else if(variable == "puJetId_simple") value = object->puJetId_simple;
940 else if(variable == "puJetId_cutbased") value = object->puJetId_cutbased;
941 else if(variable == "puJetId_tight_full") value = object->puJetId_tight_full;
942 else if(variable == "puJetId_tight_simple") value = object->puJetId_tight_simple;
943 else if(variable == "puJetId_tight_cutbased") value = object->puJetId_tight_cutbased;
944 else if(variable == "puJetId_medium_full") value = object->puJetId_medium_full;
945 else if(variable == "puJetId_medium_simple") value = object->puJetId_medium_simple;
946 else if(variable == "puJetId_medium_cutbased") value = object->puJetId_medium_cutbased;
947 else if(variable == "puJetId_loose_full") value = object->puJetId_loose_full;
948 else if(variable == "puJetId_loose_simple") value = object->puJetId_loose_simple;
949 else if(variable == "puJetId_loose_cutbased") value = object->puJetId_loose_cutbased;
950
951
952 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
953
954 value = applyFunction(function, value);
955
956 return value;
957 }
958
959
960
961 double
962 OSUAnalysis::valueLookup (const BNmuon* object, string variable, string function){
963
964 double value = 0.0;
965 if(variable == "energy") value = object->energy;
966 else if(variable == "et") value = object->et;
967 else if(variable == "pt") value = object->pt;
968 else if(variable == "px") value = object->px;
969 else if(variable == "py") value = object->py;
970 else if(variable == "pz") value = object->pz;
971 else if(variable == "phi") value = object->phi;
972 else if(variable == "eta") value = object->eta;
973 else if(variable == "theta") value = object->theta;
974 else if(variable == "trackIso") value = object->trackIso;
975 else if(variable == "ecalIso") value = object->ecalIso;
976 else if(variable == "hcalIso") value = object->hcalIso;
977 else if(variable == "caloIso") value = object->caloIso;
978 else if(variable == "trackIsDR03") value = object->trackIsoDR03;
979 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
980 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
981 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
982 else if(variable == "trackVetoIsoDR03") value = object->trackVetoIsoDR03;
983 else if(variable == "ecalVetoIsoDR03") value = object->ecalVetoIsoDR03;
984 else if(variable == "hcalVetoIsoDR03") value = object->hcalVetoIsoDR03;
985 else if(variable == "caloVetoIsoDR03") value = object->caloVetoIsoDR03;
986 else if(variable == "trackIsoDR05") value = object->trackIsoDR05;
987 else if(variable == "ecalIsoDR05") value = object->ecalIsoDR05;
988 else if(variable == "hcalIsoDR05") value = object->hcalIsoDR05;
989 else if(variable == "caloIsoDR05") value = object->caloIsoDR05;
990 else if(variable == "trackVetoIsoDR05") value = object->trackVetoIsoDR05;
991 else if(variable == "ecalVetoIsoDR05") value = object->ecalVetoIsoDR05;
992 else if(variable == "hcalVetoIsoDR05") value = object->hcalVetoIsoDR05;
993 else if(variable == "caloVetoIsoDR05") value = object->caloVetoIsoDR05;
994 else if(variable == "hcalE") value = object->hcalE;
995 else if(variable == "ecalE") value = object->ecalE;
996 else if(variable == "genET") value = object->genET;
997 else if(variable == "genPT") value = object->genPT;
998 else if(variable == "genPhi") value = object->genPhi;
999 else if(variable == "genEta") value = object->genEta;
1000 else if(variable == "genMotherET") value = object->genMotherET;
1001 else if(variable == "genMotherPT") value = object->genMotherPT;
1002 else if(variable == "genMotherPhi") value = object->genMotherPhi;
1003 else if(variable == "genMotherEta") value = object->genMotherEta;
1004 else if(variable == "vx") value = object->vx;
1005 else if(variable == "vy") value = object->vy;
1006 else if(variable == "vz") value = object->vz;
1007 else if(variable == "tkNormChi2") value = object->tkNormChi2;
1008 else if(variable == "tkPT") value = object->tkPT;
1009 else if(variable == "tkEta") value = object->tkEta;
1010 else if(variable == "tkPhi") value = object->tkPhi;
1011 else if(variable == "tkDZ") value = object->tkDZ;
1012 else if(variable == "tkD0") value = object->tkD0;
1013 else if(variable == "tkD0bs") value = object->tkD0bs;
1014 else if(variable == "tkD0err") value = object->tkD0err;
1015 else if(variable == "samNormChi2") value = object->samNormChi2;
1016 else if(variable == "samPT") value = object->samPT;
1017 else if(variable == "samEta") value = object->samEta;
1018 else if(variable == "samPhi") value = object->samPhi;
1019 else if(variable == "samDZ") value = object->samDZ;
1020 else if(variable == "samD0") value = object->samD0;
1021 else if(variable == "samD0bs") value = object->samD0bs;
1022 else if(variable == "samD0err") value = object->samD0err;
1023 else if(variable == "comNormChi2") value = object->comNormChi2;
1024 else if(variable == "comPT") value = object->comPT;
1025 else if(variable == "comEta") value = object->comEta;
1026 else if(variable == "comPhi") value = object->comPhi;
1027 else if(variable == "comDZ") value = object->comDZ;
1028 else if(variable == "comD0") value = object->comD0;
1029 else if(variable == "comD0bs") value = object->comD0bs;
1030 else if(variable == "comD0err") value = object->comD0err;
1031 else if(variable == "isolationR03emVetoEt") value = object->isolationR03emVetoEt;
1032 else if(variable == "isolationR03hadVetoEt") value = object->isolationR03hadVetoEt;
1033 else if(variable == "normalizedChi2") value = object->normalizedChi2;
1034 else if(variable == "dVzPVz") value = object->dVzPVz;
1035 else if(variable == "dB") value = object->dB;
1036 else if(variable == "ptErr") value = object->ptErr;
1037 else if(variable == "innerTrackNormChi2") value = object->innerTrackNormChi2;
1038 else if(variable == "correctedD0") value = object->correctedD0;
1039 else if(variable == "correctedD0Vertex") value = object->correctedD0Vertex;
1040 else if(variable == "correctedDZ") value = object->correctedDZ;
1041 else if(variable == "particleIso") value = object->particleIso;
1042 else if(variable == "chargedHadronIso") value = object->chargedHadronIso;
1043 else if(variable == "neutralHadronIso") value = object->neutralHadronIso;
1044 else if(variable == "photonIso") value = object->photonIso;
1045 else if(variable == "puChargedHadronIso") value = object->puChargedHadronIso;
1046 else if(variable == "chargedHadronIsoDR03") value = object->chargedHadronIsoDR03;
1047 else if(variable == "neutralHadronIsoDR03") value = object->neutralHadronIsoDR03;
1048 else if(variable == "photonIsoDR03") value = object->photonIsoDR03;
1049 else if(variable == "puChargedHadronIsoDR03") value = object->puChargedHadronIsoDR03;
1050 else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
1051 else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1052 else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1053 else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1054 else if(variable == "rhoPrime") value = object->rhoPrime;
1055 else if(variable == "AEffDr03") value = object->AEffDr03;
1056 else if(variable == "AEffDr04") value = object->AEffDr04;
1057 else if(variable == "pfIsoR03SumChargedHadronPt") value = object->pfIsoR03SumChargedHadronPt;
1058 else if(variable == "pfIsoR03SumNeutralHadronEt") value = object->pfIsoR03SumNeutralHadronEt;
1059 else if(variable == "pfIsoR03SumPhotonEt") value = object->pfIsoR03SumPhotonEt;
1060 else if(variable == "pfIsoR03SumPUPt") value = object->pfIsoR03SumPUPt;
1061 else if(variable == "pfIsoR04SumChargedHadronPt") value = object->pfIsoR04SumChargedHadronPt;
1062 else if(variable == "pfIsoR04SumNeutralHadronEt") value = object->pfIsoR04SumNeutralHadronEt;
1063 else if(variable == "pfIsoR04SumPhotonEt") value = object->pfIsoR04SumPhotonEt;
1064 else if(variable == "pfIsoR04SumPUPt") value = object->pfIsoR04SumPUPt;
1065 else if(variable == "IP") value = object->IP;
1066 else if(variable == "IPError") value = object->IPError;
1067 else if(variable == "timeAtIpInOut") value = object->timeAtIpInOut;
1068 else if(variable == "timeAtIpInOutErr") value = object->timeAtIpInOutErr;
1069 else if(variable == "timeAtIpOutIn") value = object->timeAtIpOutIn;
1070 else if(variable == "timeAtIpOutInErr") value = object->timeAtIpOutInErr;
1071 else if(variable == "ecal_time") value = object->ecal_time;
1072 else if(variable == "hcal_time") value = object->hcal_time;
1073 else if(variable == "ecal_timeError") value = object->ecal_timeError;
1074 else if(variable == "hcal_timeError") value = object->hcal_timeError;
1075 else if(variable == "energy_ecal") value = object->energy_ecal;
1076 else if(variable == "energy_hcal") value = object->energy_hcal;
1077 else if(variable == "e3x3_ecal") value = object->e3x3_ecal;
1078 else if(variable == "e3x3_hcal") value = object->e3x3_hcal;
1079 else if(variable == "energyMax_ecal") value = object->energyMax_ecal;
1080 else if(variable == "energyMax_hcal") value = object->energyMax_hcal;
1081 else if(variable == "charge") value = object->charge;
1082 else if(variable == "IDGMPTight") value = object->IDGMPTight;
1083 else if(variable == "tkNumValidHits") value = object->tkNumValidHits;
1084 else if(variable == "tkCharge") value = object->tkCharge;
1085 else if(variable == "samNumValidHits") value = object->samNumValidHits;
1086 else if(variable == "samCharge") value = object->samCharge;
1087 else if(variable == "comNumValidHits") value = object->comNumValidHits;
1088 else if(variable == "comCharge") value = object->comCharge;
1089 else if(variable == "genId") value = object->genId;
1090 else if(variable == "genCharge") value = object->genCharge;
1091 else if(variable == "genNumberOfMothers") value = object->genNumberOfMothers;
1092 else if(variable == "genMotherId") value = object->genMotherId;
1093 else if(variable == "genMotherCharge") value = object->genMotherCharge;
1094 else if(variable == "genMother0Id") value = object->genMother0Id;
1095 else if(variable == "genMother1Id") value = object->genMother1Id;
1096 else if(variable == "genGrandMother00Id") value = object->genGrandMother00Id;
1097 else if(variable == "genGrandMother01Id") value = object->genGrandMother01Id;
1098 else if(variable == "genGrandMother10Id") value = object->genGrandMother10Id;
1099 else if(variable == "genGrandMother11Id") value = object->genGrandMother11Id;
1100 else if(variable == "isPFMuon") value = object->isPFMuon;
1101 else if(variable == "isGoodMuon_1StationTight") value = object->isGoodMuon_1StationTight;
1102 else if(variable == "isGlobalMuon") value = object->isGlobalMuon;
1103 else if(variable == "isTrackerMuon") value = object->isTrackerMuon;
1104 else if(variable == "isStandAloneMuon") value = object->isStandAloneMuon;
1105 else if(variable == "isGlobalMuonPromptTight") value = object->isGlobalMuonPromptTight;
1106 else if(variable == "numberOfValidMuonHits") value = object->numberOfValidMuonHits;
1107 else if(variable == "numberOfValidTrackerHits") value = object->numberOfValidTrackerHits;
1108 else if(variable == "numberOfLayersWithMeasurement") value = object->numberOfLayersWithMeasurement;
1109 else if(variable == "pixelLayersWithMeasurement") value = object->pixelLayersWithMeasurement;
1110 else if(variable == "numberOfMatches") value = object->numberOfMatches;
1111 else if(variable == "numberOfValidTrackerHitsInnerTrack") value = object->numberOfValidTrackerHitsInnerTrack;
1112 else if(variable == "numberOfValidPixelHits") value = object->numberOfValidPixelHits;
1113 else if(variable == "numberOfMatchedStations") value = object->numberOfMatchedStations;
1114 else if(variable == "time_ndof") value = object->time_ndof;
1115
1116 //user-defined variables
1117 else if(variable == "correctedD0VertexErr") value = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1118 else if(variable == "correctedD0VertexSig") value = object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1119 else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1120 else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1121 else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1122 else if(variable == "metMT") {
1123 const BNmet *met = chosenMET ();
1124 if (met)
1125 {
1126 TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1127 p2 (met->px, met->py, 0.0, met->pt);
1128
1129 value = (p1 + p2).Mt ();
1130 }
1131 else
1132 value = -999;
1133 }
1134
1135
1136
1137 else if(variable == "correctedD0VertexInEBPlus"){
1138 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1139 else value = -999;
1140 }
1141 else if(variable == "correctedD0VertexOutEBPlus"){
1142 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1143 else value = -999;
1144 }
1145 else if(variable == "correctedD0VertexEEPlus"){
1146 if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex;
1147 else value = -999;
1148 }
1149
1150 else if(variable == "correctedD0BeamspotInEBPlus"){
1151 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1152 else value = -999;
1153 }
1154 else if(variable == "correctedD0BeamspotOutEBPlus"){
1155 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1156 else value = -999;
1157 }
1158 else if(variable == "correctedD0BeamspotEEPlus"){
1159 if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0;
1160 else value = -999;
1161 }
1162
1163 else if(variable == "correctedD0VertexInEBMinus"){
1164 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1165 else value = -999;
1166 }
1167 else if(variable == "correctedD0VertexOutEBMinus"){
1168 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1169 else value = -999;
1170 }
1171 else if(variable == "correctedD0VertexEEMinus"){
1172 if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex;
1173 else value = -999;
1174 }
1175
1176 else if(variable == "correctedD0BeamspotInEBMinus"){
1177 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1178 else value = -999;
1179 }
1180 else if(variable == "correctedD0BeamspotOutEBMinus"){
1181 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1182 else value = -999;
1183 }
1184 else if(variable == "correctedD0BeamspotEEMinus"){
1185 if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0;
1186 else value = -999;
1187 }
1188
1189
1190 else if(variable == "correctedD0VertexInEBPositiveCharge"){
1191 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1192 else value = -999;
1193 }
1194 else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1195 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1196 else value = -999;
1197 }
1198 else if(variable == "correctedD0VertexEEPositiveCharge"){
1199 if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex;
1200 else value = -999;
1201 }
1202
1203 else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1204 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1205 else value = -999;
1206 }
1207 else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1208 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1209 else value = -999;
1210 }
1211 else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1212 if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0;
1213 else value = -999;
1214 }
1215
1216 else if(variable == "correctedD0VertexInEBNegativeCharge"){
1217 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1218 else value = -999;
1219 }
1220 else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1221 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1222 else value = -999;
1223 }
1224 else if(variable == "correctedD0VertexEENegativeCharge"){
1225 if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex;
1226 else value = -999;
1227 }
1228
1229 else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1230 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1231 else value = -999;
1232 }
1233 else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1234 if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1235 else value = -999;
1236 }
1237 else if(variable == "correctedD0BeamspotEENegativeCharge"){
1238 if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0;
1239 else value = -999;
1240 }
1241
1242
1243 else if(variable == "tightID") {
1244 value = object->isGlobalMuon > 0 \
1245 && object->isPFMuon > 0 \
1246 && object->normalizedChi2 < 10 \
1247 && object->numberOfValidMuonHits > 0 \
1248 && object->numberOfMatchedStations > 1 \
1249 && fabs(object->correctedD0Vertex) < 0.2 \
1250 && fabs(object->correctedDZ) < 0.5 \
1251 && object->numberOfValidPixelHits > 0 \
1252 && object->numberOfLayersWithMeasurement > 5;
1253 }
1254 else if(variable == "tightIDdisplaced"){
1255 value = object->isGlobalMuon > 0 \
1256 && object->normalizedChi2 < 10 \
1257 && object->numberOfValidMuonHits > 0 \
1258 && object->numberOfMatchedStations > 1 \
1259 && object->numberOfValidPixelHits > 0 \
1260 && object->numberOfLayersWithMeasurement > 5;
1261 }
1262
1263 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1264
1265 else if(variable == "genMatchedPdgId"){
1266 int index = getGenMatchedParticleIndex(object);
1267 if(index == -1) value = 0;
1268 else value = mcparticles->at(index).id;
1269 }
1270
1271 else if(variable == "genMatchedId"){
1272 int index = getGenMatchedParticleIndex(object);
1273 if(index == -1) value = 0;
1274 else value = getPdgIdBinValue(mcparticles->at(index).id);
1275 }
1276 else if(variable == "genMatchedMotherId"){
1277 int index = getGenMatchedParticleIndex(object);
1278 if(index == -1) value = 0;
1279 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1280 }
1281 else if(variable == "genMatchedMotherIdReverse"){
1282 int index = getGenMatchedParticleIndex(object);
1283 if(index == -1) value = 23;
1284 else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1285 }
1286 else if(variable == "genMatchedGrandmotherId"){
1287 int index = getGenMatchedParticleIndex(object);
1288 if(index == -1) value = 0;
1289 else if(fabs(mcparticles->at(index).motherId) == 15){
1290 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1291 if(motherIndex == -1) value = 0;
1292 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1293 }
1294 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1295 }
1296
1297
1298
1299 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1300
1301 value = applyFunction(function, value);
1302
1303 return value;
1304 }
1305
1306
1307 double
1308 OSUAnalysis::valueLookup (const BNelectron* object, string variable, string function){
1309
1310 double value = 0.0;
1311 if(variable == "energy") value = object->energy;
1312 else if(variable == "et") value = object->et;
1313 else if(variable == "gsfEt") value = object->gsfEt;
1314 else if(variable == "pt") value = object->pt;
1315 else if(variable == "px") value = object->px;
1316 else if(variable == "py") value = object->py;
1317 else if(variable == "pz") value = object->pz;
1318 else if(variable == "phi") value = object->phi;
1319 else if(variable == "eta") value = object->eta;
1320 else if(variable == "theta") value = object->theta;
1321 else if(variable == "pIn") value = object->pIn;
1322 else if(variable == "pOut") value = object->pOut;
1323 else if(variable == "EscOverPin") value = object->EscOverPin;
1324 else if(variable == "EseedOverPout") value = object->EseedOverPout;
1325 else if(variable == "hadOverEm") value = object->hadOverEm;
1326 else if(variable == "trackIso") value = object->trackIso;
1327 else if(variable == "ecalIso") value = object->ecalIso;
1328 else if(variable == "hcalIso") value = object->hcalIso;
1329 else if(variable == "caloIso") value = object->caloIso;
1330 else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
1331 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
1332 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
1333 else if(variable == "hcalIsoDR03depth1") value = object->hcalIsoDR03depth1;
1334 else if(variable == "hcalIsoDR03depth2") value = object->hcalIsoDR03depth2;
1335 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
1336 else if(variable == "trackIsoDR04") value = object->trackIsoDR04;
1337 else if(variable == "ecalIsoDR04") value = object->ecalIsoDR04;
1338 else if(variable == "hcalIsoDR04") value = object->hcalIsoDR04;
1339 else if(variable == "hcalIsoDR04depth1") value = object->hcalIsoDR04depth1;
1340 else if(variable == "hcalIsoDR04depth2") value = object->hcalIsoDR04depth2;
1341 else if(variable == "caloIsoDR04") value = object->caloIsoDR04;
1342 else if(variable == "fbrem") value = object->fbrem;
1343 else if(variable == "absInvEMinusInvPin") value = object->absInvEMinusInvPin;
1344 else if(variable == "delPhiIn") value = object->delPhiIn;
1345 else if(variable == "delEtaIn") value = object->delEtaIn;
1346 else if(variable == "genET") value = object->genET;
1347 else if(variable == "genPT") value = object->genPT;
1348 else if(variable == "genPhi") value = object->genPhi;
1349 else if(variable == "genEta") value = object->genEta;
1350 else if(variable == "genMotherET") value = object->genMotherET;
1351 else if(variable == "genMotherPT") value = object->genMotherPT;
1352 else if(variable == "genMotherPhi") value = object->genMotherPhi;
1353 else if(variable == "genMotherEta") value = object->genMotherEta;
1354 else if(variable == "vx") value = object->vx;
1355 else if(variable == "vy") value = object->vy;
1356 else if(variable == "vz") value = object->vz;
1357 else if(variable == "scEnergy") value = object->scEnergy;
1358 else if(variable == "scRawEnergy") value = object->scRawEnergy;
1359 else if(variable == "scSigmaEtaEta") value = object->scSigmaEtaEta;
1360 else if(variable == "scSigmaIEtaIEta") value = object->scSigmaIEtaIEta;
1361 else if(variable == "scE1x5") value = object->scE1x5;
1362 else if(variable == "scE2x5Max") value = object->scE2x5Max;
1363 else if(variable == "scE5x5") value = object->scE5x5;
1364 else if(variable == "scEt") value = object->scEt;
1365 else if(variable == "scEta") value = object->scEta;
1366 else if(variable == "scPhi") value = object->scPhi;
1367 else if(variable == "scZ") value = object->scZ;
1368 else if(variable == "tkNormChi2") value = object->tkNormChi2;
1369 else if(variable == "tkPT") value = object->tkPT;
1370 else if(variable == "tkEta") value = object->tkEta;
1371 else if(variable == "tkPhi") value = object->tkPhi;
1372 else if(variable == "tkDZ") value = object->tkDZ;
1373 else if(variable == "tkD0") value = object->tkD0;
1374 else if(variable == "tkD0bs") value = object->tkD0bs;
1375 else if(variable == "tkD0err") value = object->tkD0err;
1376 else if(variable == "mva") value = object->mva;
1377 else if(variable == "mvaTrigV0") value = object->mvaTrigV0;
1378 else if(variable == "mvaNonTrigV0") value = object->mvaNonTrigV0;
1379 else if(variable == "dist") value = object->dist;
1380 else if(variable == "dcot") value = object->dcot;
1381 else if(variable == "convradius") value = object->convradius;
1382 else if(variable == "convPointX") value = object->convPointX;
1383 else if(variable == "convPointY") value = object->convPointY;
1384 else if(variable == "convPointZ") value = object->convPointZ;
1385 else if(variable == "eMax") value = object->eMax;
1386 else if(variable == "eLeft") value = object->eLeft;
1387 else if(variable == "eRight") value = object->eRight;
1388 else if(variable == "eTop") value = object->eTop;
1389 else if(variable == "eBottom") value = object->eBottom;
1390 else if(variable == "e3x3") value = object->e3x3;
1391 else if(variable == "swissCross") value = object->swissCross;
1392 else if(variable == "seedEnergy") value = object->seedEnergy;
1393 else if(variable == "seedTime") value = object->seedTime;
1394 else if(variable == "swissCrossNoI85") value = object->swissCrossNoI85;
1395 else if(variable == "swissCrossI85") value = object->swissCrossI85;
1396 else if(variable == "E2overE9NoI85") value = object->E2overE9NoI85;
1397 else if(variable == "E2overE9I85") value = object->E2overE9I85;
1398 else if(variable == "correctedD0") value = object->correctedD0;
1399 else if(variable == "correctedD0Vertex") value = object->correctedD0Vertex;
1400 else if(variable == "correctedDZ") value = object->correctedDZ;
1401 else if(variable == "particleIso") value = object->particleIso;
1402 else if(variable == "chargedHadronIso") value = object->chargedHadronIso;
1403 else if(variable == "neutralHadronIso") value = object->neutralHadronIso;
1404 else if(variable == "photonIso") value = object->photonIso;
1405 else if(variable == "puChargedHadronIso") value = object->puChargedHadronIso;
1406 else if(variable == "chargedHadronIsoDR03") value = object->chargedHadronIsoDR03;
1407 else if(variable == "neutralHadronIsoDR03") value = object->neutralHadronIsoDR03;
1408 else if(variable == "photonIsoDR03") value = object->photonIsoDR03;
1409 else if(variable == "puChargedHadronIsoDR03") value = object->puChargedHadronIsoDR03;
1410 else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
1411 else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1412 else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1413 else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1414 else if(variable == "rhoPrime") value = object->rhoPrime;
1415 else if(variable == "AEffDr03") value = object->AEffDr03;
1416 else if(variable == "AEffDr04") value = object->AEffDr04;
1417 else if(variable == "IP") value = object->IP;
1418 else if(variable == "IPError") value = object->IPError;
1419 else if(variable == "charge") value = object->charge;
1420 else if(variable == "classification") value = object->classification;
1421 else if(variable == "genId") value = object->genId;
1422 else if(variable == "genCharge") value = object->genCharge;
1423 else if(variable == "genNumberOfMothers") value = object->genNumberOfMothers;
1424 else if(variable == "genMotherId") value = object->genMotherId;
1425 else if(variable == "genMotherCharge") value = object->genMotherCharge;
1426 else if(variable == "genMother0Id") value = object->genMother0Id;
1427 else if(variable == "genMother1Id") value = object->genMother1Id;
1428 else if(variable == "genGrandMother00Id") value = object->genGrandMother00Id;
1429 else if(variable == "genGrandMother01Id") value = object->genGrandMother01Id;
1430 else if(variable == "genGrandMother10Id") value = object->genGrandMother10Id;
1431 else if(variable == "genGrandMother11Id") value = object->genGrandMother11Id;
1432 else if(variable == "numClusters") value = object->numClusters;
1433 else if(variable == "tkNumValidHits") value = object->tkNumValidHits;
1434 else if(variable == "tkCharge") value = object->tkCharge;
1435 else if(variable == "gsfCharge") value = object->gsfCharge;
1436 else if(variable == "isEB") value = object->isEB;
1437 else if(variable == "isEE") value = object->isEE;
1438 else if(variable == "isGap") value = object->isGap;
1439 else if(variable == "isEBEEGap") value = object->isEBEEGap;
1440 else if(variable == "isEBGap") value = object->isEBGap;
1441 else if(variable == "isEEGap") value = object->isEEGap;
1442 else if(variable == "isEcalDriven") value = object->isEcalDriven;
1443 else if(variable == "isTrackerDriven") value = object->isTrackerDriven;
1444 else if(variable == "numberOfLostHits") value = object->numberOfLostHits;
1445 else if(variable == "numberOfExpectedInnerHits") value = object->numberOfExpectedInnerHits;
1446 else if(variable == "numberOfValidPixelHits") value = object->numberOfValidPixelHits;
1447 else if(variable == "numberOfValidPixelBarrelHits") value = object->numberOfValidPixelBarrelHits;
1448 else if(variable == "numberOfValidPixelEndcapHits") value = object->numberOfValidPixelEndcapHits;
1449 else if(variable == "isHEEP") value = object->isHEEP;
1450 else if(variable == "isHEEPnoEt") value = object->isHEEPnoEt;
1451 else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
1452 else if(variable == "eidRobustHighEnergy") value = object->eidRobustHighEnergy;
1453 else if(variable == "eidRobustLoose") value = object->eidRobustLoose;
1454 else if(variable == "eidRobustTight") value = object->eidRobustTight;
1455 else if(variable == "eidLoose") value = object->eidLoose;
1456 else if(variable == "eidTight") value = object->eidTight;
1457 else if(variable == "eidVeryLooseMC") value = object->eidVeryLooseMC;
1458 else if(variable == "eidLooseMC") value = object->eidLooseMC;
1459 else if(variable == "eidMediumMC") value = object->eidMediumMC;
1460 else if(variable == "eidTightMC") value = object->eidTightMC;
1461 else if(variable == "eidSuperTightMC") value = object->eidSuperTightMC;
1462 else if(variable == "eidHyperTight1MC") value = object->eidHyperTight1MC;
1463 else if(variable == "eidHyperTight2MC") value = object->eidHyperTight2MC;
1464 else if(variable == "eidHyperTight3MC") value = object->eidHyperTight3MC;
1465 else if(variable == "eidHyperTight4MC") value = object->eidHyperTight4MC;
1466 else if(variable == "passConvVeto") value = object->passConvVeto;
1467
1468 //user-defined variables
1469 else if(variable == "correctedD0VertexErr") value = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1470 else if(variable == "correctedD0VertexSig") value = object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1471 else if(variable == "detIso") value = (object->trackIso) / object->pt;
1472 else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1473 else if(variable == "metMT") {
1474 const BNmet *met = chosenMET ();
1475 if (met)
1476 {
1477 TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1478 p2 (met->px, met->py, 0.0, met->pt);
1479
1480 value = (p1 + p2).Mt ();
1481 }
1482 else
1483 value = -999;
1484 }
1485
1486 else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
1487 if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
1488 else value = -999;
1489 }
1490 else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){
1491 if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex;
1492 else value = -999;
1493 }
1494
1495 else if(variable == "correctedD0VertexInEBPlus"){
1496 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1497 else value = -999;
1498 }
1499 else if(variable == "correctedD0VertexOutEBPlus"){
1500 if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1501 else value = -999;
1502 }
1503 else if(variable == "correctedD0VertexEEPlus"){
1504 if(object->isEE && object->eta > 0) value = object->correctedD0Vertex;
1505 else value = -999;
1506 }
1507
1508 else if(variable == "correctedD0BeamspotInEBPlus"){
1509 if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1510 else value = -999;
1511 }
1512 else if(variable == "correctedD0BeamspotOutEBPlus"){
1513 if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1514 else value = -999;
1515 }
1516 else if(variable == "correctedD0BeamspotEEPlus"){
1517 if(object->isEE && object->eta > 0) value = object->correctedD0;
1518 else value = -999;
1519 }
1520
1521 else if(variable == "correctedD0VertexInEBMinus"){
1522 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1523 else value = -999;
1524 }
1525 else if(variable == "correctedD0VertexOutEBMinus"){
1526 if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1527 else value = -999;
1528 }
1529 else if(variable == "correctedD0VertexEEMinus"){
1530 if(object->isEE && object->eta < 0) value = object->correctedD0Vertex;
1531 else value = -999;
1532 }
1533
1534 else if(variable == "correctedD0BeamspotInEBMinus"){
1535 if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1536 else value = -999;
1537 }
1538 else if(variable == "correctedD0BeamspotOutEBMinus"){
1539 if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1540 else value = -999;
1541 }
1542 else if(variable == "correctedD0BeamspotEEMinus"){
1543 if(object->isEE && object->eta < 0) value = object->correctedD0;
1544 else value = -999;
1545 }
1546
1547 else if(variable == "tightID"){
1548 if (object->isEB)
1549 {
1550 value = fabs(object->delEtaIn) < 0.004 \
1551 && fabs (object->delPhiIn) < 0.03 \
1552 && object->scSigmaIEtaIEta < 0.01 \
1553 && object->hadOverEm < 0.12 \
1554 && fabs (object->correctedD0Vertex) < 0.02 \
1555 && fabs (object->correctedDZ) < 0.1 \
1556 && object->absInvEMinusInvPin < 0.05 \
1557 && object->passConvVeto;
1558 }
1559 else
1560 {
1561 value = fabs(object->delEtaIn) < 0.005 \
1562 && fabs (object->delPhiIn) < 0.02 \
1563 && object->scSigmaIEtaIEta < 0.03 \
1564 && object->hadOverEm < 0.10 \
1565 && fabs (object->correctedD0Vertex) < 0.02 \
1566 && fabs (object->correctedDZ) < 0.1 \
1567 && object->absInvEMinusInvPin < 0.05 \
1568 && object->passConvVeto;
1569 }
1570 }
1571
1572 else if(variable == "correctedD0VertexInEBPositiveCharge"){
1573 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1574 else value = -999;
1575 }
1576 else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1577 if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1578 else value = -999;
1579 }
1580 else if(variable == "correctedD0VertexEEPositiveCharge"){
1581 if(object->isEE && object->charge > 0) value = object->correctedD0Vertex;
1582 else value = -999;
1583 }
1584
1585 else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1586 if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1587 else value = -999;
1588 }
1589 else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1590 if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1591 else value = -999;
1592 }
1593 else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1594 if(object->isEE && object->charge > 0) value = object->correctedD0;
1595 else value = -999;
1596 }
1597
1598 else if(variable == "correctedD0VertexInEBNegativeCharge"){
1599 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1600 else value = -999;
1601 }
1602 else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1603 if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1604 else value = -999;
1605 }
1606 else if(variable == "correctedD0VertexEENegativeCharge"){
1607 if(object->isEE && object->charge < 0) value = object->correctedD0Vertex;
1608 else value = -999;
1609 }
1610
1611 else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1612 if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1613 else value = -999;
1614 }
1615 else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1616 if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1617 else value = -999;
1618 }
1619 else if(variable == "correctedD0BeamspotEENegativeCharge"){
1620 if(object->isEE && object->charge < 0) value = object->correctedD0;
1621 else value = -999;
1622 }
1623
1624
1625 else if(variable == "tightIDdisplaced"){
1626 if (object->isEB)
1627 {
1628 value = fabs(object->delEtaIn) < 0.004 \
1629 && fabs (object->delPhiIn) < 0.03 \
1630 && object->scSigmaIEtaIEta < 0.01 \
1631 && object->hadOverEm < 0.12 \
1632 && object->absInvEMinusInvPin < 0.05;
1633 }
1634 else
1635 {
1636 value = fabs (object->delEtaIn) < 0.005 \
1637 && fabs (object->delPhiIn) < 0.02 \
1638 && object->scSigmaIEtaIEta < 0.03 \
1639 && object->hadOverEm < 0.10 \
1640 && object->absInvEMinusInvPin < 0.05;
1641 }
1642 }
1643
1644
1645 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1646
1647 else if(variable == "genMatchedPdgId"){
1648 int index = getGenMatchedParticleIndex(object);
1649 if(index == -1) value = 0;
1650 else value = mcparticles->at(index).id;
1651 }
1652
1653
1654 else if(variable == "genMatchedId"){
1655 int index = getGenMatchedParticleIndex(object);
1656 if(index == -1) value = 0;
1657 else value = getPdgIdBinValue(mcparticles->at(index).id);
1658 }
1659 else if(variable == "genMatchedMotherId"){
1660 int index = getGenMatchedParticleIndex(object);
1661 if(index == -1) value = 0;
1662 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1663 }
1664 else if(variable == "genMatchedMotherIdReverse"){
1665 int index = getGenMatchedParticleIndex(object);
1666 if(index == -1) value = 23;
1667 else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1668 }
1669 else if(variable == "genMatchedGrandmotherId"){
1670 int index = getGenMatchedParticleIndex(object);
1671 if(index == -1) value = 0;
1672 else if(fabs(mcparticles->at(index).motherId) == 15){
1673 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1674 if(motherIndex == -1) value = 0;
1675 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1676 }
1677 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1678 }
1679
1680
1681
1682 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1683
1684 value = applyFunction(function, value);
1685
1686 return value;
1687 }
1688
1689
1690 double
1691 OSUAnalysis::valueLookup (const BNevent* object, string variable, string function){
1692
1693 double value = 0.0;
1694
1695 if(variable == "weight") value = object->weight;
1696 else if(variable == "pthat") value = object->pthat;
1697 else if(variable == "qScale") value = object->qScale;
1698 else if(variable == "alphaQCD") value = object->alphaQCD;
1699 else if(variable == "alphaQED") value = object->alphaQED;
1700 else if(variable == "scalePDF") value = object->scalePDF;
1701 else if(variable == "x1") value = object->x1;
1702 else if(variable == "x2") value = object->x2;
1703 else if(variable == "xPDF1") value = object->xPDF1;
1704 else if(variable == "xPDF2") value = object->xPDF2;
1705 else if(variable == "BSx") value = object->BSx;
1706 else if(variable == "BSy") value = object->BSy;
1707 else if(variable == "BSz") value = object->BSz;
1708 else if(variable == "PVx") value = object->PVx;
1709 else if(variable == "PVy") value = object->PVy;
1710 else if(variable == "PVz") value = object->PVz;
1711 else if(variable == "bField") value = object->bField;
1712 else if(variable == "instLumi") value = object->instLumi;
1713 else if(variable == "bxLumi") value = object->bxLumi;
1714 else if(variable == "FilterOutScrapingFraction") value = object->FilterOutScrapingFraction;
1715 else if(variable == "sumNVtx") value = object->sumNVtx;
1716 else if(variable == "sumTrueNVtx") value = object->sumTrueNVtx;
1717 else if(variable == "nm1_true") value = object->nm1_true;
1718 else if(variable == "n0_true") value = object->n0_true;
1719 else if(variable == "np1_true") value = object->np1_true;
1720 else if(variable == "numTruePV") value = object->numTruePV;
1721 else if(variable == "Q2ScaleUpWgt") value = object->Q2ScaleUpWgt;
1722 else if(variable == "Q2ScaleDownWgt") value = object->Q2ScaleDownWgt;
1723 else if(variable == "rho_kt6PFJets") value = object->rho_kt6PFJets;
1724 else if(variable == "rho_kt6PFJetsCentralChargedPileUp") value = object->rho_kt6PFJetsCentralChargedPileUp;
1725 else if(variable == "rho_kt6PFJetsCentralNeutral") value = object->rho_kt6PFJetsCentralNeutral;
1726 else if(variable == "rho_kt6PFJetsCentralNeutralTight") value = object->rho_kt6PFJetsCentralNeutralTight;
1727 else if(variable == "run") value = object->run;
1728 else if(variable == "lumi") value = object->lumi;
1729 else if(variable == "sample") value = object->sample;
1730 else if(variable == "numPV") value = object->numPV;
1731 else if(variable == "W0decay") value = object->W0decay;
1732 else if(variable == "W1decay") value = object->W1decay;
1733 else if(variable == "Z0decay") value = object->Z0decay;
1734 else if(variable == "Z1decay") value = object->Z1decay;
1735 else if(variable == "H0decay") value = object->H0decay;
1736 else if(variable == "H1decay") value = object->H1decay;
1737 else if(variable == "hcalnoiseLoose") value = object->hcalnoiseLoose;
1738 else if(variable == "hcalnoiseTight") value = object->hcalnoiseTight;
1739 else if(variable == "GoodVertex") value = object->GoodVertex;
1740 else if(variable == "FilterOutScraping") value = object->FilterOutScraping;
1741 else if(variable == "HBHENoiseFilter") value = object->HBHENoiseFilter;
1742 else if(variable == "CSCLooseHaloId") value = object->CSCLooseHaloId;
1743 else if(variable == "CSCTightHaloId") value = object->CSCTightHaloId;
1744 else if(variable == "EcalLooseHaloId") value = object->EcalLooseHaloId;
1745 else if(variable == "EcalTightHaloId") value = object->EcalTightHaloId;
1746 else if(variable == "HcalLooseHaloId") value = object->HcalLooseHaloId;
1747 else if(variable == "HcalTightHaloId") value = object->HcalTightHaloId;
1748 else if(variable == "GlobalLooseHaloId") value = object->GlobalLooseHaloId;
1749 else if(variable == "GlobalTightHaloId") value = object->GlobalTightHaloId;
1750 else if(variable == "LooseId") value = object->LooseId;
1751 else if(variable == "TightId") value = object->TightId;
1752 else if(variable == "numGenPV") value = object->numGenPV;
1753 else if(variable == "nm1") value = object->nm1;
1754 else if(variable == "n0") value = object->n0;
1755 else if(variable == "np1") value = object->np1;
1756 else if(variable == "id1") value = object->id1;
1757 else if(variable == "id2") value = object->id2;
1758 else if(variable == "evt") value = object->evt;
1759 else if(variable == "puScaleFactor"){
1760 if(datasetType_ != "data")
1761 value = puWeight_->at (events->at (0).numTruePV);
1762 else
1763 value = 1.0;
1764 }
1765 else if(variable == "muonScaleFactor"){
1766 if(datasetType_ != "data")
1767 // value = muonSFWeight_->at (chosenMuon ()->eta);
1768 value = 1.0;
1769 else
1770 value = 1.0;
1771 }
1772 else if(variable == "electronScaleFactor"){
1773 if(datasetType_ != "data")
1774 // value = electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
1775 value = 1.0;
1776 else
1777 value = 1.0;
1778 }
1779
1780 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1781
1782 value = applyFunction(function, value);
1783
1784 return value;
1785 }
1786
1787 double
1788 OSUAnalysis::valueLookup (const BNtau* object, string variable, string function){
1789
1790 double value = 0.0;
1791
1792 if(variable == "px") value = object->px;
1793 else if(variable == "py") value = object->py;
1794 else if(variable == "pz") value = object->pz;
1795 else if(variable == "energy") value = object->energy;
1796 else if(variable == "et") value = object->et;
1797 else if(variable == "pt") value = object->pt;
1798 else if(variable == "eta") value = object->eta;
1799 else if(variable == "phi") value = object->phi;
1800 else if(variable == "emFraction") value = object->emFraction;
1801 else if(variable == "leadingTrackPt") value = object->leadingTrackPt;
1802 else if(variable == "leadingTrackIpVtdxy") value = object->leadingTrackIpVtdxy;
1803 else if(variable == "leadingTrackIpVtdz") value = object->leadingTrackIpVtdz;
1804 else if(variable == "leadingTrackIpVtdxyError") value = object->leadingTrackIpVtdxyError;
1805 else if(variable == "leadingTrackIpVtdzError") value = object->leadingTrackIpVtdzError;
1806 else if(variable == "leadingTrackVx") value = object->leadingTrackVx;
1807 else if(variable == "leadingTrackVy") value = object->leadingTrackVy;
1808 else if(variable == "leadingTrackVz") value = object->leadingTrackVz;
1809 else if(variable == "leadingTrackValidHits") value = object->leadingTrackValidHits;
1810 else if(variable == "leadingTrackNormChiSqrd") value = object->leadingTrackNormChiSqrd;
1811 else if(variable == "numProngs") value = object->numProngs;
1812 else if(variable == "numSignalGammas") value = object->numSignalGammas;
1813 else if(variable == "numSignalNeutrals") value = object->numSignalNeutrals;
1814 else if(variable == "numSignalPiZeros") value = object->numSignalPiZeros;
1815 else if(variable == "decayMode") value = object->decayMode;
1816 else if(variable == "charge") value = object->charge;
1817 else if(variable == "inTheCracks") value = object->inTheCracks;
1818 else if(variable == "HPSagainstElectronLoose") value = object->HPSagainstElectronLoose;
1819 else if(variable == "HPSagainstElectronMVA") value = object->HPSagainstElectronMVA;
1820 else if(variable == "HPSagainstElectronMedium") value = object->HPSagainstElectronMedium;
1821 else if(variable == "HPSagainstElectronTight") value = object->HPSagainstElectronTight;
1822 else if(variable == "HPSagainstMuonLoose") value = object->HPSagainstMuonLoose;
1823 else if(variable == "HPSagainstMuonMedium") value = object->HPSagainstMuonMedium;
1824 else if(variable == "HPSagainstMuonTight") value = object->HPSagainstMuonTight;
1825 else if(variable == "HPSbyLooseCombinedIsolationDeltaBetaCorr") value = object->HPSbyLooseCombinedIsolationDeltaBetaCorr;
1826 else if(variable == "HPSbyMediumCombinedIsolationDeltaBetaCorr") value = object->HPSbyMediumCombinedIsolationDeltaBetaCorr;
1827 else if(variable == "HPSbyTightCombinedIsolationDeltaBetaCorr") value = object->HPSbyTightCombinedIsolationDeltaBetaCorr;
1828 else if(variable == "HPSbyVLooseCombinedIsolationDeltaBetaCorr") value = object->HPSbyVLooseCombinedIsolationDeltaBetaCorr;
1829 else if(variable == "HPSdecayModeFinding") value = object->HPSdecayModeFinding;
1830 else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1831
1832
1833
1834 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1835
1836 else if(variable == "genMatchedPdgId"){
1837 int index = getGenMatchedParticleIndex(object);
1838 if(index == -1) value = 0;
1839 else value = mcparticles->at(index).id;
1840 }
1841
1842 else if(variable == "genMatchedId"){
1843 int index = getGenMatchedParticleIndex(object);
1844 if(index == -1) value = 0;
1845 else value = getPdgIdBinValue(mcparticles->at(index).id);
1846 }
1847 else if(variable == "genMatchedMotherId"){
1848 int index = getGenMatchedParticleIndex(object);
1849 if(index == -1) value = 0;
1850 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1851 }
1852 else if(variable == "genMatchedMotherIdReverse"){
1853 int index = getGenMatchedParticleIndex(object);
1854 if(index == -1) value = 23;
1855 else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1856 }
1857 else if(variable == "genMatchedGrandmotherId"){
1858 int index = getGenMatchedParticleIndex(object);
1859 if(index == -1) value = 0;
1860 else if(fabs(mcparticles->at(index).motherId) == 15){
1861 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1862 if(motherIndex == -1) value = 0;
1863 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1864 }
1865 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1866 }
1867
1868
1869 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1870
1871 value = applyFunction(function, value);
1872
1873 return value;
1874 }
1875
1876 double
1877 OSUAnalysis::valueLookup (const BNmet* object, string variable, string function){
1878
1879 double value = 0.0;
1880
1881 if(variable == "et") value = object->et;
1882 else if(variable == "pt") value = object->pt;
1883 else if(variable == "px") value = object->px;
1884 else if(variable == "py") value = object->py;
1885 else if(variable == "phi") value = object->phi;
1886 else if(variable == "Upt") value = object->Upt;
1887 else if(variable == "Uphi") value = object->Uphi;
1888 else if(variable == "NeutralEMFraction") value = object->NeutralEMFraction;
1889 else if(variable == "NeutralHadEtFraction") value = object->NeutralHadEtFraction;
1890 else if(variable == "ChargedEMEtFraction") value = object->ChargedEMEtFraction;
1891 else if(variable == "ChargedHadEtFraction") value = object->ChargedHadEtFraction;
1892 else if(variable == "MuonEtFraction") value = object->MuonEtFraction;
1893 else if(variable == "Type6EtFraction") value = object->Type6EtFraction;
1894 else if(variable == "Type7EtFraction") value = object->Type7EtFraction;
1895 else if(variable == "genPT") value = object->genPT;
1896 else if(variable == "genPhi") value = object->genPhi;
1897 else if(variable == "muonCorEx") value = object->muonCorEx;
1898 else if(variable == "muonCorEy") value = object->muonCorEy;
1899 else if(variable == "jet20CorEx") value = object->jet20CorEx;
1900 else if(variable == "jet20CorEy") value = object->jet20CorEy;
1901 else if(variable == "jet1CorEx") value = object->jet1CorEx;
1902 else if(variable == "jet1CorEy") value = object->jet1CorEy;
1903 else if(variable == "sumET") value = object->sumET;
1904 else if(variable == "corSumET") value = object->corSumET;
1905 else if(variable == "mEtSig") value = object->mEtSig;
1906 else if(variable == "metSignificance") value = object->metSignificance;
1907 else if(variable == "significance") value = object->significance;
1908 else if(variable == "sigmaX2") value = object->sigmaX2;
1909 else if(variable == "sigmaY2") value = object->sigmaY2;
1910 else if(variable == "sigmaXY") value = object->sigmaXY;
1911 else if(variable == "sigmaYX") value = object->sigmaYX;
1912 else if(variable == "maxEtInEmTowers") value = object->maxEtInEmTowers;
1913 else if(variable == "emEtFraction") value = object->emEtFraction;
1914 else if(variable == "emEtInEB") value = object->emEtInEB;
1915 else if(variable == "emEtInEE") value = object->emEtInEE;
1916 else if(variable == "emEtInHF") value = object->emEtInHF;
1917 else if(variable == "maxEtInHadTowers") value = object->maxEtInHadTowers;
1918 else if(variable == "hadEtFraction") value = object->hadEtFraction;
1919 else if(variable == "hadEtInHB") value = object->hadEtInHB;
1920 else if(variable == "hadEtInHE") value = object->hadEtInHE;
1921 else if(variable == "hadEtInHF") value = object->hadEtInHF;
1922 else if(variable == "hadEtInHO") value = object->hadEtInHO;
1923 else if(variable == "UDeltaPx") value = object->UDeltaPx;
1924 else if(variable == "UDeltaPy") value = object->UDeltaPy;
1925 else if(variable == "UDeltaP") value = object->UDeltaP;
1926 else if(variable == "Uscale") value = object->Uscale;
1927 else if(variable == "type2corPx") value = object->type2corPx;
1928 else if(variable == "type2corPy") value = object->type2corPy;
1929 else if(variable == "T2pt") value = object->T2pt;
1930 else if(variable == "T2px") value = object->T2px;
1931 else if(variable == "T2py") value = object->T2py;
1932 else if(variable == "T2phi") value = object->T2phi;
1933 else if(variable == "T2sumET") value = object->T2sumET;
1934 else if(variable == "pfT1jet1pt") value = object->pfT1jet1pt;
1935 else if(variable == "pfT1jet1phi") value = object->pfT1jet1phi;
1936 else if(variable == "pfT1jet6pt") value = object->pfT1jet6pt;
1937 else if(variable == "pfT1jet6phi") value = object->pfT1jet6phi;
1938 else if(variable == "pfT1jet10pt") value = object->pfT1jet10pt;
1939 else if(variable == "pfT1jet10phi") value = object->pfT1jet10phi;
1940
1941 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1942
1943 value = applyFunction(function, value);
1944
1945 return value;
1946 }
1947
1948 double
1949 OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){
1950
1951 double value = 0.0;
1952 double pMag = sqrt(object->pt * object->pt +
1953 object->pz * object->pz);
1954
1955 if(variable == "pt") value = object->pt;
1956 else if(variable == "px") value = object->px;
1957 else if(variable == "py") value = object->py;
1958 else if(variable == "pz") value = object->pz;
1959 else if(variable == "phi") value = object->phi;
1960 else if(variable == "eta") value = object->eta;
1961 else if(variable == "theta") value = object->theta;
1962 else if(variable == "normChi2") value = object->normChi2;
1963 else if(variable == "dZ") value = object->dZ;
1964 else if(variable == "d0") value = object->d0;
1965 else if(variable == "d0err") value = object->d0err;
1966 else if(variable == "vx") value = object->vx;
1967 else if(variable == "vy") value = object->vy;
1968 else if(variable == "vz") value = object->vz;
1969 else if(variable == "charge") value = object->charge;
1970 else if(variable == "numValidHits") value = object->numValidHits;
1971 else if(variable == "isHighPurity") value = object->isHighPurity;
1972
1973
1974 //additional BNs info for disappTrks
1975 else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
1976 else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
1977 else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
1978 else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
1979 else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
1980 else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
1981 else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
1982 else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
1983 else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
1984 else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
1985
1986
1987 //user defined variables
1988 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;
1989 else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
1990 else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
1991 else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
1992 else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
1993 else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
1994 else if(variable == "ptErrorByPt") value = (object->ptError/object->pt);
1995 else if(variable == "ptError") value = object->ptError;
1996 else if(variable == "ptRes") value = getTrkPtRes(object);
1997
1998 else if (variable == "d0wrtPV"){
1999 double vx = object->vx - chosenVertex ()->x,
2000 vy = object->vy - chosenVertex ()->y,
2001 px = object->px,
2002 py = object->py,
2003 pt = object->pt;
2004 value = (-vx * py + vy * px) / pt;
2005 }
2006 else if (variable == "dZwrtPV"){
2007 double vx = object->vx - chosenVertex ()->x,
2008 vy = object->vy - chosenVertex ()->y,
2009 vz = object->vz - chosenVertex ()->z,
2010 px = object->px,
2011 py = object->py,
2012 pz = object->pz,
2013 pt = object->pt;
2014 value = vz - (vx * px + vy * py)/pt * (pz/pt);
2015 }
2016
2017
2018
2019 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2020
2021 else if(variable == "genMatchedPdgId"){
2022 int index = getGenMatchedParticleIndex(object);
2023 if(index == -1) value = 0;
2024 else value = mcparticles->at(index).id;
2025 }
2026
2027 else if(variable == "genMatchedId"){
2028 int index = getGenMatchedParticleIndex(object);
2029 if(index == -1) value = 0;
2030 else value = getPdgIdBinValue(mcparticles->at(index).id);
2031 }
2032 else if(variable == "genMatchedMotherId"){
2033 int index = getGenMatchedParticleIndex(object);
2034 if(index == -1) value = 0;
2035 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2036 }
2037 else if(variable == "genMatchedMotherIdReverse"){
2038 int index = getGenMatchedParticleIndex(object);
2039 if(index == -1) value = 23;
2040 else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2041 }
2042 else if(variable == "genMatchedGrandmotherId"){
2043 int index = getGenMatchedParticleIndex(object);
2044 if(index == -1) value = 0;
2045 else if(fabs(mcparticles->at(index).motherId) == 15){
2046 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2047 if(motherIndex == -1) value = 0;
2048 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2049 }
2050 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2051 }
2052
2053
2054
2055 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2056
2057 value = applyFunction(function, value);
2058
2059 return value;
2060 }
2061
2062 double
2063 OSUAnalysis::valueLookup (const BNgenjet* object, string variable, string function){
2064
2065 double value = 0.0;
2066
2067 if(variable == "pt") value = object->pt;
2068 else if(variable == "eta") value = object->eta;
2069 else if(variable == "phi") value = object->phi;
2070 else if(variable == "px") value = object->px;
2071 else if(variable == "py") value = object->py;
2072 else if(variable == "pz") value = object->pz;
2073 else if(variable == "et") value = object->et;
2074 else if(variable == "energy") value = object->energy;
2075 else if(variable == "mass") value = object->mass;
2076 else if(variable == "emEnergy") value = object->emEnergy;
2077 else if(variable == "hadEnergy") value = object->hadEnergy;
2078 else if(variable == "invisibleEnergy") value = object->invisibleEnergy;
2079 else if(variable == "auxiliaryEnergy") value = object->auxiliaryEnergy;
2080 else if(variable == "charge") value = object->charge;
2081
2082
2083 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2084
2085 value = applyFunction(function, value);
2086
2087 return value;
2088 }
2089
2090 double
2091 OSUAnalysis::valueLookup (const BNmcparticle* object, string variable, string function){
2092
2093 double value = 0.0;
2094
2095 if(variable == "energy") value = object->energy;
2096 else if(variable == "et") value = object->et;
2097 else if(variable == "pt") value = object->pt;
2098 else if(variable == "px") value = object->px;
2099 else if(variable == "py") value = object->py;
2100 else if(variable == "pz") value = object->pz;
2101 else if(variable == "phi") value = object->phi;
2102 else if(variable == "eta") value = object->eta;
2103 else if(variable == "theta") value = object->theta;
2104 else if(variable == "mass") value = object->mass;
2105 else if(variable == "vx") value = object->vx;
2106 else if(variable == "vy") value = object->vy;
2107 else if(variable == "vz") value = object->vz;
2108 else if(variable == "motherET") value = object->motherET;
2109 else if(variable == "motherPT") value = object->motherPT;
2110 else if(variable == "motherPhi") value = object->motherPhi;
2111 else if(variable == "motherEta") value = object->motherEta;
2112 else if(variable == "mother0ET") value = object->mother0ET;
2113 else if(variable == "mother0PT") value = object->mother0PT;
2114 else if(variable == "mother0Phi") value = object->mother0Phi;
2115 else if(variable == "mother0Eta") value = object->mother0Eta;
2116 else if(variable == "mother1ET") value = object->mother1ET;
2117 else if(variable == "mother1PT") value = object->mother1PT;
2118 else if(variable == "mother1Phi") value = object->mother1Phi;
2119 else if(variable == "mother1Eta") value = object->mother1Eta;
2120 else if(variable == "daughter0ET") value = object->daughter0ET;
2121 else if(variable == "daughter0PT") value = object->daughter0PT;
2122 else if(variable == "daughter0Phi") value = object->daughter0Phi;
2123 else if(variable == "daughter0Eta") value = object->daughter0Eta;
2124 else if(variable == "daughter1ET") value = object->daughter1ET;
2125 else if(variable == "daughter1PT") value = object->daughter1PT;
2126 else if(variable == "daughter1Phi") value = object->daughter1Phi;
2127 else if(variable == "daughter1Eta") value = object->daughter1Eta;
2128 else if(variable == "grandMotherET") value = object->grandMotherET;
2129 else if(variable == "grandMotherPT") value = object->grandMotherPT;
2130 else if(variable == "grandMotherPhi") value = object->grandMotherPhi;
2131 else if(variable == "grandMotherEta") value = object->grandMotherEta;
2132 else if(variable == "grandMother00ET") value = object->grandMother00ET;
2133 else if(variable == "grandMother00PT") value = object->grandMother00PT;
2134 else if(variable == "grandMother00Phi") value = object->grandMother00Phi;
2135 else if(variable == "grandMother00Eta") value = object->grandMother00Eta;
2136 else if(variable == "grandMother01ET") value = object->grandMother01ET;
2137 else if(variable == "grandMother01PT") value = object->grandMother01PT;
2138 else if(variable == "grandMother01Phi") value = object->grandMother01Phi;
2139 else if(variable == "grandMother01Eta") value = object->grandMother01Eta;
2140 else if(variable == "grandMother10ET") value = object->grandMother10ET;
2141 else if(variable == "grandMother10PT") value = object->grandMother10PT;
2142 else if(variable == "grandMother10Phi") value = object->grandMother10Phi;
2143 else if(variable == "grandMother10Eta") value = object->grandMother10Eta;
2144 else if(variable == "grandMother11ET") value = object->grandMother11ET;
2145 else if(variable == "grandMother11PT") value = object->grandMother11PT;
2146 else if(variable == "grandMother11Phi") value = object->grandMother11Phi;
2147 else if(variable == "grandMother11Eta") value = object->grandMother11Eta;
2148 else if(variable == "charge") value = object->charge;
2149 else if(variable == "id") value = object->id;
2150 else if(variable == "status") value = object->status;
2151 else if(variable == "motherId") value = object->motherId;
2152 else if(variable == "motherCharge") value = object->motherCharge;
2153 else if(variable == "mother0Id") value = object->mother0Id;
2154 else if(variable == "mother0Status") value = object->mother0Status;
2155 else if(variable == "mother0Charge") value = object->mother0Charge;
2156 else if(variable == "mother1Id") value = object->mother1Id;
2157 else if(variable == "mother1Status") value = object->mother1Status;
2158 else if(variable == "mother1Charge") value = object->mother1Charge;
2159 else if(variable == "daughter0Id") value = object->daughter0Id;
2160 else if(variable == "daughter0Status") value = object->daughter0Status;
2161 else if(variable == "daughter0Charge") value = object->daughter0Charge;
2162 else if(variable == "daughter1Id") value = object->daughter1Id;
2163 else if(variable == "daughter1Status") value = object->daughter1Status;
2164 else if(variable == "daughter1Charge") value = object->daughter1Charge;
2165 else if(variable == "grandMotherId") value = object->grandMotherId;
2166 else if(variable == "grandMotherCharge") value = object->grandMotherCharge;
2167 else if(variable == "grandMother00Id") value = object->grandMother00Id;
2168 else if(variable == "grandMother00Status") value = object->grandMother00Status;
2169 else if(variable == "grandMother00Charge") value = object->grandMother00Charge;
2170 else if(variable == "grandMother01Id") value = object->grandMother01Id;
2171 else if(variable == "grandMother01Status") value = object->grandMother01Status;
2172 else if(variable == "grandMother01Charge") value = object->grandMother01Charge;
2173 else if(variable == "grandMother10Id") value = object->grandMother10Id;
2174 else if(variable == "grandMother10Status") value = object->grandMother10Status;
2175 else if(variable == "grandMother10Charge") value = object->grandMother10Charge;
2176 else if(variable == "grandMother11Id") value = object->grandMother11Id;
2177 else if(variable == "grandMother11Status") value = object->grandMother11Status;
2178 else if(variable == "grandMother11Charge") value = object->grandMother11Charge;
2179
2180 //user-defined variables
2181 else if (variable == "d0"){
2182 double vx = object->vx - chosenVertex ()->x,
2183 vy = object->vy - chosenVertex ()->y,
2184 px = object->px,
2185 py = object->py,
2186 pt = object->pt;
2187 value = (-vx * py + vy * px) / pt;
2188 }
2189 else if (variable == "dz"){
2190 double vx = object->vx - chosenVertex ()->x,
2191 vy = object->vy - chosenVertex ()->y,
2192 vz = object->vz - chosenVertex ()->z,
2193 px = object->px,
2194 py = object->py,
2195 pz = object->pz,
2196 pt = object->pt;
2197 value = vz - (vx * px + vy * py)/pt * (pz/pt);
2198 }
2199 else if(variable == "v0"){
2200 value = sqrt(object->vx*object->vx + object->vy*object->vy);
2201 }
2202 else if(variable == "deltaV0"){
2203 value = sqrt((object->vx-chosenVertex ()->x)*(object->vx-chosenVertex ()->x) + (object->vy-chosenVertex ()->y)*(object->vy-chosenVertex ()->y));
2204 }
2205 else if (variable == "deltaVx"){
2206 value = object->vx - chosenVertex ()->x;
2207 }
2208 else if (variable == "deltaVy"){
2209 value = object->vy - chosenVertex ()->y;
2210 }
2211 else if (variable == "deltaVz"){
2212 value = object->vz - chosenVertex ()->z;
2213 }
2214
2215
2216 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2217
2218 value = applyFunction(function, value);
2219
2220 return value;
2221 }
2222
2223 double
2224 OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable, string function){
2225
2226 double value = 0.0;
2227
2228 if(variable == "x") value = object->x;
2229 else if(variable == "xError") value = object->xError;
2230 else if(variable == "y") value = object->y;
2231 else if(variable == "yError") value = object->yError;
2232 else if(variable == "z") value = object->z;
2233 else if(variable == "zError") value = object->zError;
2234 else if(variable == "rho") value = object->rho;
2235 else if(variable == "normalizedChi2") value = object->normalizedChi2;
2236 else if(variable == "ndof") value = object->ndof;
2237 else if(variable == "isFake") value = object->isFake;
2238 else if(variable == "isValid") value = object->isValid;
2239 else if(variable == "tracksSize") value = object->tracksSize;
2240 else if(variable == "isGood") value = object->isGood;
2241
2242
2243 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2244
2245 value = applyFunction(function, value);
2246
2247 return value;
2248 }
2249
2250 double
2251 OSUAnalysis::valueLookup (const BNbxlumi* object, string variable, string function){
2252
2253 double value = 0.0;
2254
2255 if(variable == "bx_B1_now") value = object->bx_B1_now;
2256 else if(variable == "bx_B2_now") value = object->bx_B2_now;
2257 else if(variable == "bx_LUMI_now") value = object->bx_LUMI_now;
2258
2259
2260 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2261
2262 value = applyFunction(function, value);
2263
2264 return value;
2265 }
2266
2267 double
2268 OSUAnalysis::valueLookup (const BNphoton* object, string variable, string function){
2269
2270 double value = 0.0;
2271
2272 if(variable == "energy") value = object->energy;
2273 else if(variable == "et") value = object->et;
2274 else if(variable == "pt") value = object->pt;
2275 else if(variable == "px") value = object->px;
2276 else if(variable == "py") value = object->py;
2277 else if(variable == "pz") value = object->pz;
2278 else if(variable == "phi") value = object->phi;
2279 else if(variable == "eta") value = object->eta;
2280 else if(variable == "theta") value = object->theta;
2281 else if(variable == "trackIso") value = object->trackIso;
2282 else if(variable == "ecalIso") value = object->ecalIso;
2283 else if(variable == "hcalIso") value = object->hcalIso;
2284 else if(variable == "caloIso") value = object->caloIso;
2285 else if(variable == "trackIsoHollowConeDR03") value = object->trackIsoHollowConeDR03;
2286 else if(variable == "trackIsoSolidConeDR03") value = object->trackIsoSolidConeDR03;
2287 else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
2288 else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
2289 else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
2290 else if(variable == "trackIsoHollowConeDR04") value = object->trackIsoHollowConeDR04;
2291 else if(variable == "trackIsoSolidConeDR04") value = object->trackIsoSolidConeDR04;
2292 else if(variable == "ecalIsoDR04") value = object->ecalIsoDR04;
2293 else if(variable == "hcalIsoDR04") value = object->hcalIsoDR04;
2294 else if(variable == "caloIsoDR04") value = object->caloIsoDR04;
2295 else if(variable == "hadOverEm") value = object->hadOverEm;
2296 else if(variable == "sigmaEtaEta") value = object->sigmaEtaEta;
2297 else if(variable == "sigmaIetaIeta") value = object->sigmaIetaIeta;
2298 else if(variable == "r9") value = object->r9;
2299 else if(variable == "scEnergy") value = object->scEnergy;
2300 else if(variable == "scRawEnergy") value = object->scRawEnergy;
2301 else if(variable == "scSeedEnergy") value = object->scSeedEnergy;
2302 else if(variable == "scEta") value = object->scEta;
2303 else if(variable == "scPhi") value = object->scPhi;
2304 else if(variable == "scZ") value = object->scZ;
2305 else if(variable == "genET") value = object->genET;
2306 else if(variable == "genPT") value = object->genPT;
2307 else if(variable == "genPhi") value = object->genPhi;
2308 else if(variable == "genEta") value = object->genEta;
2309 else if(variable == "genMotherET") value = object->genMotherET;
2310 else if(variable == "genMotherPT") value = object->genMotherPT;
2311 else if(variable == "genMotherPhi") value = object->genMotherPhi;
2312 else if(variable == "genMotherEta") value = object->genMotherEta;
2313 else if(variable == "eMax") value = object->eMax;
2314 else if(variable == "eLeft") value = object->eLeft;
2315 else if(variable == "eRight") value = object->eRight;
2316 else if(variable == "eTop") value = object->eTop;
2317 else if(variable == "eBottom") value = object->eBottom;
2318 else if(variable == "e3x3") value = object->e3x3;
2319 else if(variable == "swissCross") value = object->swissCross;
2320 else if(variable == "seedEnergy") value = object->seedEnergy;
2321 else if(variable == "seedTime") value = object->seedTime;
2322 else if(variable == "swissCrossNoI85") value = object->swissCrossNoI85;
2323 else if(variable == "swissCrossI85") value = object->swissCrossI85;
2324 else if(variable == "E2overE9NoI85") value = object->E2overE9NoI85;
2325 else if(variable == "E2overE9I85") value = object->E2overE9I85;
2326 else if(variable == "IDTight") value = object->IDTight;
2327 else if(variable == "IDLoose") value = object->IDLoose;
2328 else if(variable == "IDLooseEM") value = object->IDLooseEM;
2329 else if(variable == "genId") value = object->genId;
2330 else if(variable == "genCharge") value = object->genCharge;
2331 else if(variable == "genMotherId") value = object->genMotherId;
2332 else if(variable == "genMotherCharge") value = object->genMotherCharge;
2333 else if(variable == "isEB") value = object->isEB;
2334 else if(variable == "isEE") value = object->isEE;
2335 else if(variable == "isGap") value = object->isGap;
2336 else if(variable == "isEBEEGap") value = object->isEBEEGap;
2337 else if(variable == "isEBGap") value = object->isEBGap;
2338 else if(variable == "isEEGap") value = object->isEEGap;
2339 else if(variable == "hasPixelSeed") value = object->hasPixelSeed;
2340 else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
2341
2342
2343 else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2344
2345 else if(variable == "genMatchedPdgId"){
2346 int index = getGenMatchedParticleIndex(object);
2347 if(index == -1) value = 0;
2348 else value = mcparticles->at(index).id;
2349 }
2350 else if(variable == "genMatchedId"){
2351 int index = getGenMatchedParticleIndex(object);
2352 if(index == -1) value = 0;
2353 else value = getPdgIdBinValue(mcparticles->at(index).id);
2354 }
2355 else if(variable == "genMatchedMotherId"){
2356 int index = getGenMatchedParticleIndex(object);
2357 if(index == -1) value = 0;
2358 else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2359 }
2360 else if(variable == "genMatchedMotherIdReverse"){
2361 int index = getGenMatchedParticleIndex(object);
2362 if(index == -1) value = 23;
2363 else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2364 }
2365 else if(variable == "genMatchedGrandmotherId"){
2366 int index = getGenMatchedParticleIndex(object);
2367 if(index == -1) value = 0;
2368 else if(fabs(mcparticles->at(index).motherId) == 15){
2369 int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2370 if(motherIndex == -1) value = 0;
2371 else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2372 }
2373 else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2374 }
2375
2376
2377 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2378
2379 value = applyFunction(function, value);
2380
2381 return value;
2382 }
2383
2384 double
2385 OSUAnalysis::valueLookup (const BNsupercluster* object, string variable, string function){
2386
2387 double value = 0.0;
2388
2389 if(variable == "energy") value = object->energy;
2390 else if(variable == "et") value = object->et;
2391 else if(variable == "ex") value = object->ex;
2392 else if(variable == "ey") value = object->ey;
2393 else if(variable == "ez") value = object->ez;
2394 else if(variable == "phi") value = object->phi;
2395 else if(variable == "eta") value = object->eta;
2396 else if(variable == "theta") value = object->theta;
2397
2398
2399 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2400
2401 value = applyFunction(function, value);
2402
2403 return value;
2404 }
2405
2406
2407 double
2408 OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function){
2409
2410 double value = 0.0;
2411
2412 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2413 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2414 else if(variable == "invMass"){
2415 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2416 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2417 value = (fourVector1 + fourVector2).M();
2418 }
2419 else if(variable == "pt"){
2420 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2421 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2422 value = (fourVector1 + fourVector2).Et();
2423 }
2424 else if(variable == "threeDAngle")
2425 {
2426 TVector3 threeVector1(object1->px, object1->py, object1->pz);
2427 TVector3 threeVector2(object2->px, object2->py, object2->pz);
2428 value = (threeVector1.Angle(threeVector2));
2429 }
2430 else if(variable == "alpha")
2431 {
2432 static const double pi = 3.1415926535897932384626433832795028841971693993751058;
2433 TVector3 threeVector1(object1->px, object1->py, object1->pz);
2434 TVector3 threeVector2(object2->px, object2->py, object2->pz);
2435 value = (pi-threeVector1.Angle(threeVector2));
2436 }
2437 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2438 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2439 else if(variable == "d0Sign"){
2440 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2441 if(d0Sign < 0) value = -0.5;
2442 else if (d0Sign > 0) value = 0.5;
2443 else value = -999;
2444 }
2445 else if(variable == "chargeProduct"){
2446 value = object1->charge*object2->charge;
2447 }
2448 else if(variable == "muon1CorrectedD0Vertex"){
2449 value = object1->correctedD0Vertex;
2450 }
2451 else if(variable == "muon2CorrectedD0Vertex"){
2452 value = object2->correctedD0Vertex;
2453 }
2454 else if(variable == "muon1timeAtIpInOut"){
2455 value = object1->timeAtIpInOut;
2456 }
2457 else if(variable == "muon2timeAtIpInOut"){
2458 value = object2->timeAtIpInOut;
2459 }
2460 else if(variable == "muon1correctedD0")
2461 {
2462 value = object1->correctedD0;
2463 }
2464 else if(variable == "muon2correctedD0")
2465 {
2466 value = object2->correctedD0;
2467 }
2468
2469 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2470
2471 value = applyFunction(function, value);
2472
2473 return value;
2474 }
2475
2476 double
2477 OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function){
2478
2479 double value = 0.0;
2480
2481 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2482 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2483 else if(variable == "invMass"){
2484 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2485 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2486 value = (fourVector1 + fourVector2).M();
2487 }
2488 else if(variable == "pt"){
2489 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2490 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2491 value = (fourVector1 + fourVector2).Et();
2492 }
2493 else if(variable == "threeDAngle")
2494 {
2495 TVector3 threeVector1(object1->px, object1->py, object1->pz);
2496 TVector3 threeVector2(object2->px, object2->py, object2->pz);
2497 value = (threeVector1.Angle(threeVector2));
2498 }
2499 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2500 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2501 else if(variable == "d0Sign"){
2502 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2503 if(d0Sign < 0) value = -0.5;
2504 else if (d0Sign > 0) value = 0.5;
2505 else value = -999;
2506 }
2507 else if(variable == "chargeProduct"){
2508 value = object1->charge*object2->charge;
2509 }
2510 else if(variable == "electron1CorrectedD0Vertex"){
2511 value = object1->correctedD0Vertex;
2512 }
2513 else if(variable == "electron2CorrectedD0Vertex"){
2514 value = object2->correctedD0Vertex;
2515 }
2516 else if(variable == "electron1CorrectedD0"){
2517 value = object1->correctedD0;
2518 }
2519 else if(variable == "electron2CorrectedD0"){
2520 value = object2->correctedD0;
2521 }
2522
2523 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2524
2525 value = applyFunction(function, value);
2526
2527 return value;
2528 }
2529
2530 double
2531 OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function){
2532
2533 double value = 0.0;
2534
2535 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2536 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2537 else if(variable == "invMass"){
2538 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2539 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2540 value = (fourVector1 + fourVector2).M();
2541 }
2542 else if(variable == "pt"){
2543 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2544 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2545 value = (fourVector1 + fourVector2).Et();
2546 }
2547 else if(variable == "threeDAngle")
2548 {
2549 TVector3 threeVector1(object1->px, object1->py, object1->pz);
2550 TVector3 threeVector2(object2->px, object2->py, object2->pz);
2551 value = (threeVector1.Angle(threeVector2));
2552 }
2553 else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2554 else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2555 else if(variable == "d0Sign"){
2556 double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2557 if(d0Sign < 0) value = -0.5;
2558 else if (d0Sign > 0) value = 0.5;
2559 else value = -999;
2560 }
2561 else if(variable == "chargeProduct"){
2562 value = object1->charge*object2->charge;
2563 }
2564 else if(variable == "electronCorrectedD0Vertex"){
2565 value = object1->correctedD0Vertex;
2566 }
2567 else if(variable == "muonCorrectedD0Vertex"){
2568 value = object2->correctedD0Vertex;
2569 }
2570 else if(variable == "electronCorrectedD0"){
2571 value = object1->correctedD0;
2572 }
2573 else if(variable == "muonCorrectedD0"){
2574 value = object2->correctedD0;
2575 }
2576 else if(variable == "electronDetIso"){
2577 value = (object1->trackIso) / object1->pt;
2578 }
2579 else if(variable == "muonDetIso"){
2580 value = (object2->trackIsoDR03) / object2->pt;
2581 }
2582
2583
2584 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2585
2586 value = applyFunction(function, value);
2587
2588 return value;
2589 }
2590
2591
2592 double
2593 OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrack* object2, string variable, string function){
2594 double electronMass = 0.000511;
2595 double value = 0.0;
2596 TLorentzVector fourVector1(0, 0, 0, 0);
2597 TLorentzVector fourVector2(0, 0, 0, 0);
2598 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2599 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2600 else if(variable == "invMass"){
2601 fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, electronMass);
2602 fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, electronMass );
2603
2604 value = (fourVector1 + fourVector2).M();
2605 }
2606
2607 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2608 value = applyFunction(function, value);
2609 return value;
2610 }
2611
2612 double
2613 OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrack* object2, string variable, string function){
2614 double pionMass = 0.140;
2615 double muonMass = 0.106;
2616 double value = 0.0;
2617 TLorentzVector fourVector1(0, 0, 0, 0);
2618 TLorentzVector fourVector2(0, 0, 0, 0);
2619 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2620 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2621 else if(variable == "invMass"){
2622 fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, muonMass);
2623 fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, pionMass );
2624
2625 value = (fourVector1 + fourVector2).M();
2626 }
2627
2628 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2629 value = applyFunction(function, value);
2630 return value;
2631 }
2632
2633
2634 double
2635 OSUAnalysis::valueLookup (const BNtau* object1, const BNtau* object2, string variable, string function){
2636 double value = 0.0;
2637 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2638 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2639 else if(variable == "invMass"){
2640 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2641 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2642 value = (fourVector1 + fourVector2).M();
2643 }
2644
2645 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2646 value = applyFunction(function, value);
2647 return value;
2648 }
2649
2650 double
2651 OSUAnalysis::valueLookup (const BNmuon* object1, const BNtau* object2, string variable, string function){
2652 double value = 0.0;
2653 if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2654 else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2655 else if(variable == "invMass"){
2656 TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2657 TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2658 value = (fourVector1 + fourVector2).M();
2659 }
2660
2661 else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2662 value = applyFunction(function, value);
2663 return value;
2664 }
2665
2666
2667
2668 // Calculate the number of tracks in cone of DeltaR<0.5 around track1.
2669 // Return true iff no other tracks are found in this cone.
2670 int
2671 OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
2672 for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
2673 if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.
2674 double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
2675 if(deltaRtrk < 0.5) return 0;
2676 }
2677 return 1;
2678
2679 }
2680
2681
2682 double
2683 OSUAnalysis::getTrkPtRes (const BNtrack* track1){
2684
2685 double ptTrue = getTrkPtTrue(track1, mcparticles.product());
2686 double PtRes = (track1->pt - ptTrue) / ptTrue;
2687
2688 return PtRes;
2689
2690 }
2691
2692
2693 double
2694 OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
2695 double value = -99;
2696 double genDeltaRLowest = 999;
2697
2698 for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
2699 double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
2700 if (genDeltaRtemp < genDeltaRLowest) {
2701 genDeltaRLowest = genDeltaRtemp;
2702 if (genDeltaRLowest < 0.05) { // Only consider it truth-matched if DeltaR<0.15.
2703 double ptTrue = genPart->pt;
2704 value = ptTrue;
2705 }
2706 }
2707 }
2708
2709 return value;
2710
2711 }
2712
2713 //creates a map of the dead Ecal channels in the barrel and endcap
2714 //to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
2715 //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
2716 void
2717 OSUAnalysis::WriteDeadEcal (){
2718 double etaEcal, phiEcal;
2719 ifstream DeadEcalFile(deadEcalFile_);
2720 if(!DeadEcalFile) {
2721 cout << "Error: DeadEcalFile has not been found." << endl;
2722 return;
2723 }
2724 if(DeadEcalVec.size()!= 0){
2725 cout << "Error: DeadEcalVec has a nonzero size" << endl;
2726 return;
2727 }
2728 while(!DeadEcalFile.eof())
2729 {
2730 DeadEcalFile >> etaEcal >> phiEcal;
2731 DeadEcal newChan;
2732 newChan.etaEcal = etaEcal;
2733 newChan.phiEcal = phiEcal;
2734 DeadEcalVec.push_back(newChan);
2735 }
2736 if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl;
2737 }
2738
2739 //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
2740 int
2741 OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
2742 double deltaRLowest = 999;
2743 int value = 0;
2744 if (DeadEcalVec.size() == 0) WriteDeadEcal();
2745 for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
2746 double eta = ecal->etaEcal;
2747 double phi = ecal->phiEcal;
2748 double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
2749 if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
2750 }
2751 if (deltaRLowest<0.05) {value = 1;}
2752 else {value = 0;}
2753 return value;
2754 }
2755
2756
2757 // Returns the smallest DeltaR between the object and any generated true particle in the event.
2758 template <class InputObject>
2759 double OSUAnalysis::getGenDeltaRLowest(InputObject object){
2760 double genDeltaRLowest = 999.;
2761 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2762 double deltaRtemp = deltaR(mcparticle->eta, mcparticle->phi, object->eta, object->phi);
2763 if (deltaRtemp < genDeltaRLowest) genDeltaRLowest = deltaRtemp;
2764 }
2765 return genDeltaRLowest;
2766 }
2767
2768
2769 double
2770 OSUAnalysis::applyFunction(string function, double value){
2771
2772 if(function == "abs") value = fabs(value);
2773 else if(function == "fabs") value = fabs(value);
2774 else if(function == "log10") value = log10(value);
2775 else if(function == "log") value = log10(value);
2776
2777 else if(function == "") value = value;
2778 else{std::cout << "WARNING: invalid function '" << function << "'\n";}
2779
2780 return value;
2781
2782 }
2783
2784
2785 template <class InputCollection>
2786 void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2787
2788
2789 for (uint object = 0; object != inputCollection->size(); object++){
2790
2791
2792 bool decision = true;//object passes if this cut doesn't cut on that type of object
2793
2794
2795 if(currentCut.inputCollection == inputType){
2796
2797 vector<bool> subcutDecisions;
2798 for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2799 double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2800 subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2801 }
2802 if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2803 else{
2804 bool tempDecision = true;
2805 for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2806 if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2807 tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2808 else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2809 tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
2810 }
2811 decision = tempDecision;
2812 }
2813 }
2814 individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2815
2816
2817 //set flags for objects that pass each cut AND all the previous cuts
2818 bool previousCumulativeFlag = true;
2819 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2820 if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(object)) previousCumulativeFlag = true;
2821 else{ previousCumulativeFlag = false; break;}
2822 }
2823 cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2824
2825
2826 }
2827
2828 }
2829
2830
2831 template <class InputCollection1, class InputCollection2>
2832 void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
2833 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
2834
2835
2836 bool sameObjects = false;
2837 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2838
2839
2840 int counter = 0;
2841 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2842 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2843
2844 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2845
2846
2847 bool decision = true;//object passes if this cut doesn't cut on that type of object
2848
2849 if(currentCut.inputCollection == inputType){
2850
2851 vector<bool> subcutDecisions;
2852 for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2853 double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2854 subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2855 }
2856
2857 if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2858 else{
2859 bool tempDecision = subcutDecisions.at(0);
2860 for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
2861 if(currentCut.logicalOperators.at(subcutIndex-1) == "&" || currentCut.logicalOperators.at(subcutIndex-1) == "&&")
2862 tempDecision = tempDecision && subcutDecisions.at(subcutIndex);
2863 else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
2864 tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
2865 }
2866 decision = tempDecision;
2867 }
2868 }
2869 individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2870
2871 //set flags for objects that pass each cut AND all the previous cuts
2872 bool previousCumulativeFlag = true;
2873 for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2874 if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
2875 else{ previousCumulativeFlag = false; break;}
2876 }
2877 //apply flags for the components of the composite object as well
2878 bool currentCumulativeFlag = true;
2879 if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
2880 else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
2881 else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
2882 else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
2883 cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
2884
2885 counter++;
2886 }
2887
2888 }
2889
2890
2891 }
2892
2893
2894 template <class InputCollection>
2895 void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
2896
2897 for (uint object = 0; object != inputCollection->size(); object++){
2898
2899 if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2900
2901 string currentString = parameters.inputVariables.at(0);
2902 string inputVariable = "";
2903 string function = "";
2904 if(currentString.find("(")==std::string::npos){
2905 inputVariable = currentString;// variable to cut on
2906 }
2907 else{
2908 function = currentString.substr(0,currentString.find("("));//function comes before the "("
2909 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2910 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2911 }
2912
2913 double value = valueLookup(&inputCollection->at(object), inputVariable, function);
2914 histo->Fill(value,scaleFactor);
2915
2916 if (printEventInfo_) {
2917 // Write information about event to screen, for testing purposes.
2918 cout << " Info for event: value for histogram " << histo->GetName() << ": " << value << endl;
2919 }
2920
2921 }
2922 }
2923
2924 template <class InputCollection1, class InputCollection2>
2925 void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
2926
2927 bool sameObjects = false;
2928 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2929
2930 int pairCounter = 0;
2931 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2932 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2933
2934 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2935
2936 //only take objects which have passed all cuts and pairs which have passed all cuts
2937 if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2938 if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2939 if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2940
2941 string currentString = parameters.inputVariables.at(0);
2942 string inputVariable = "";
2943 string function = "";
2944 if(currentString.find("(")==std::string::npos){
2945 inputVariable = currentString;// variable to cut on
2946 }
2947 else{
2948 function = currentString.substr(0,currentString.find("("));//function comes before the "("
2949 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2950 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2951 }
2952
2953 double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2954 histo->Fill(value,scaleFactor);
2955
2956 pairCounter++;
2957 }
2958 }
2959
2960 }
2961
2962
2963 template <class InputCollection>
2964 void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
2965
2966 for (uint object = 0; object != inputCollection->size(); object++){
2967
2968 if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2969
2970 string currentString = parameters.inputVariables.at(0);
2971 string inputVariable = "";
2972 string function = "";
2973 if(currentString.find("(")==std::string::npos){
2974 inputVariable = currentString;// variable to cut on
2975 }
2976 else{
2977 function = currentString.substr(0,currentString.find("("));//function comes before the "("
2978 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2979 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2980 }
2981 double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
2982
2983 currentString = parameters.inputVariables.at(1);
2984 inputVariable = "";
2985 function = "";
2986 if(currentString.find("(")==std::string::npos){
2987 inputVariable = currentString;// variable to cut on
2988 }
2989 else{
2990 function = currentString.substr(0,currentString.find("("));//function comes before the "("
2991 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2992 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2993 }
2994
2995 double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
2996
2997 histo->Fill(valueX,valueY,scaleFactor);
2998
2999 }
3000
3001 }
3002
3003 template <class InputCollection1, class InputCollection2>
3004 void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
3005
3006 bool sameObjects = false;
3007 if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
3008
3009 int pairCounter = 0;
3010 for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
3011 for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
3012
3013 if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
3014
3015 //only take objects which have passed all cuts and pairs which have passed all cuts
3016 if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
3017 if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
3018 if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
3019
3020 string currentString = parameters.inputVariables.at(0);
3021 string inputVariable = "";
3022 string function = "";
3023 if(currentString.find("(")==std::string::npos){
3024 inputVariable = currentString;// variable to cut on
3025 }
3026 else{
3027 function = currentString.substr(0,currentString.find("("));//function comes before the "("
3028 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3029 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3030 }
3031 double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
3032
3033 currentString = parameters.inputVariables.at(1);
3034 inputVariable = "";
3035 function = "";
3036 if(currentString.find("(")==std::string::npos){
3037 inputVariable = currentString;// variable to cut on
3038 }
3039 else{
3040 function = currentString.substr(0,currentString.find("("));//function comes before the "("
3041 inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3042 inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3043 }
3044 double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
3045
3046
3047 histo->Fill(valueX,valueY,scaleFactor);
3048
3049 pairCounter++;
3050
3051 }
3052 }
3053
3054 }
3055
3056
3057 template <class InputObject>
3058 int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
3059
3060 int bestMatchIndex = -1;
3061 double bestMatchDeltaR = 999;
3062
3063 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
3064
3065 double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
3066 if(currentDeltaR > 0.05) continue;
3067 // cout << std::setprecision(3) << std::setw(20)
3068 // << "\tcurrentParticle: eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
3069 // << std::setw(20)
3070 // << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
3071 // << std::setw(20)
3072 // << "\tdeltaR = " << currentDeltaR
3073 // << std::setprecision(1)
3074 // << std::setw(20)
3075 // << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
3076 // << std::setw(20)
3077 // << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
3078 // << std::setw(20)
3079 // << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
3080 if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
3081 bestMatchIndex = mcparticle - mcparticles->begin();
3082 bestMatchDeltaR = currentDeltaR;
3083 }
3084
3085 }
3086 // if(bestMatchDeltaR != 999) cout << "bestMatch: deltaR = " << bestMatchDeltaR << " id = " << mcparticles->at(bestMatchIndex).id << " motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
3087 // else cout << "no match found..." << endl;
3088 return bestMatchIndex;
3089
3090 }
3091
3092
3093 int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
3094
3095 int bestMatchIndex = -1;
3096 double bestMatchDeltaR = 999;
3097
3098 for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
3099
3100 if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
3101
3102 double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
3103 if(currentDeltaR > 0.05) continue;
3104
3105 if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
3106 bestMatchIndex = mcparticle - mcparticles->begin();
3107 bestMatchDeltaR = currentDeltaR;
3108 }
3109
3110 }
3111
3112 return bestMatchIndex;
3113 }
3114
3115
3116 // bin particle type
3117 // --- -------------
3118 // 0 unmatched
3119 // 1 u
3120 // 2 d
3121 // 3 s
3122 // 4 c
3123 // 5 b
3124 // 6 t
3125 // 7 e
3126 // 8 mu
3127 // 9 tau
3128 // 10 nu
3129 // 11 g
3130 // 12 gamma
3131 // 13 Z
3132 // 14 W
3133 // 15 light meson
3134 // 16 K meson
3135 // 17 D meson
3136 // 18 B meson
3137 // 19 light baryon
3138 // 20 strange baryon
3139 // 21 charm baryon
3140 // 22 bottom baryon
3141 // 23 other
3142
3143 int OSUAnalysis::getPdgIdBinValue(int pdgId){
3144
3145 int binValue = -999;
3146 int absPdgId = fabs(pdgId);
3147 if(pdgId == -1) binValue = 0;
3148 else if(absPdgId <= 6 ) binValue = absPdgId;
3149 else if(absPdgId == 11 ) binValue = 7;
3150 else if(absPdgId == 13 ) binValue = 8;
3151 else if(absPdgId == 15 ) binValue = 9;
3152 else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
3153 else if(absPdgId == 21 ) binValue = 11;
3154 else if(absPdgId == 22 ) binValue = 12;
3155 else if(absPdgId == 23 ) binValue = 13;
3156 else if(absPdgId == 24 ) binValue = 14;
3157 else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130 ) binValue = 15;
3158 else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400) ) binValue = 16;
3159 else if(absPdgId > 400 && absPdgId < 500 ) binValue = 17;
3160 else if(absPdgId > 500 && absPdgId < 600 ) binValue = 18;
3161 else if(absPdgId > 1000 && absPdgId < 3000 ) binValue = 19;
3162 else if(absPdgId > 3000 && absPdgId < 4000 ) binValue = 20;
3163 else if(absPdgId > 4000 && absPdgId < 5000 ) binValue = 21;
3164 else if(absPdgId > 5000 && absPdgId < 6000 ) binValue = 22;
3165
3166 else binValue = 23;
3167
3168 return binValue;
3169
3170 }
3171
3172 const BNprimaryvertex *
3173 OSUAnalysis::chosenVertex ()
3174 {
3175 const BNprimaryvertex *chosenVertex = 0;
3176 if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
3177 vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back().size() ? cumulativeFlags.at("primaryvertexs").back() :
3178 cumulativeFlags.at("primaryvertexs").at(cumulativeFlags.at("primaryvertexs").size() - 2);
3179 for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
3180 if(!vertexFlags.at(vertexIndex)) continue;
3181 chosenVertex = & primaryvertexs->at(vertexIndex);
3182 break;
3183 }
3184 }
3185 else {
3186 chosenVertex = &primaryvertexs->at(0);
3187 }
3188
3189 return chosenVertex;
3190 }
3191
3192 const BNmet *
3193 OSUAnalysis::chosenMET ()
3194 {
3195 const BNmet *chosenMET = 0;
3196 if(std::find(objectsToCut.begin(), objectsToCut.end(), "mets") != objectsToCut.end()) {
3197 vector<bool> metFlags = cumulativeFlags.at("mets").back().size() ? cumulativeFlags.at("mets").back() :
3198 cumulativeFlags.at("mets").at(cumulativeFlags.at("mets").size() - 2);
3199 for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
3200 if(!metFlags.at(metIndex)) continue;
3201 chosenMET = & mets->at(metIndex);
3202 break;
3203 }
3204 }
3205 else {
3206 chosenMET = &mets->at(0);
3207 }
3208
3209 return chosenMET;
3210 }
3211
3212 const BNelectron *
3213 OSUAnalysis::chosenElectron ()
3214 {
3215 const BNelectron *chosenElectron = 0;
3216 if(std::find(objectsToCut.begin(), objectsToCut.end(), "electrons") != objectsToCut.end()) {
3217 vector<bool> electronFlags = cumulativeFlags.at("electrons").back().size() ? cumulativeFlags.at("electrons").back() :
3218 cumulativeFlags.at("electrons").at(cumulativeFlags.at("electrons").size() - 2);
3219 for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
3220 if(!electronFlags.at(electronIndex)) continue;
3221 chosenElectron = & electrons->at(electronIndex);
3222 break;
3223 }
3224 }
3225 else {
3226 chosenElectron = &electrons->at(0);
3227 }
3228
3229 return chosenElectron;
3230 }
3231
3232 const BNmuon *
3233 OSUAnalysis::chosenMuon ()
3234 {
3235 const BNmuon *chosenMuon = 0;
3236 if(std::find(objectsToCut.begin(), objectsToCut.end(), "muons") != objectsToCut.end()) {
3237 vector<bool> muonFlags = cumulativeFlags.at("muons").back().size() ? cumulativeFlags.at("muons").back() :
3238 cumulativeFlags.at("muons").at(cumulativeFlags.at("muons").size() - 2);
3239 for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
3240 if(!muonFlags.at(muonIndex)) continue;
3241 chosenMuon = & muons->at(muonIndex);
3242 break;
3243 }
3244 }
3245 else {
3246 chosenMuon = &muons->at(0);
3247 }
3248
3249 return chosenMuon;
3250 }
3251
3252
3253 DEFINE_FWK_MODULE(OSUAnalysis);