ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc
(Generate patch)

Comparing UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc (file contents):
Revision 1.26 by lantonel, Mon Mar 11 13:34:26 2013 UTC vs.
Revision 1.49 by wulsin, Wed Apr 17 15:39:43 2013 UTC

# Line 11 | Line 11 | OSUAnalysis::OSUAnalysis (const edm::Par
11    tracks_ (cfg.getParameter<edm::InputTag> ("tracks")),
12    genjets_ (cfg.getParameter<edm::InputTag> ("genjets")),
13    mcparticles_ (cfg.getParameter<edm::InputTag> ("mcparticles")),
14 +  stops_ (cfg.getParameter<edm::InputTag> ("stops")),
15    primaryvertexs_ (cfg.getParameter<edm::InputTag> ("primaryvertexs")),
16    bxlumis_ (cfg.getParameter<edm::InputTag> ("bxlumis")),
17    photons_ (cfg.getParameter<edm::InputTag> ("photons")),
# Line 18 | Line 19 | OSUAnalysis::OSUAnalysis (const edm::Par
19    triggers_ (cfg.getParameter<edm::InputTag> ("triggers")),
20    puFile_ (cfg.getParameter<std::string> ("puFile")),
21    deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")),
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"))
30 >  plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
31 >  doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting")),
32 >  printEventInfo_ (cfg.getParameter<bool> ("printEventInfo")),
33 >  useTrackCaloRhoCorr_ (cfg.getParameter<bool> ("useTrackCaloRhoCorr")),
34 >  stopCTau_ (cfg.getParameter<vector<double> > ("stopCTau"))
35  
36   {
37  
38    TH1::SetDefaultSumw2 ();
39  
40    //create pile-up reweighting object, if necessary
41 <  if(datasetType_ != "data") puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
41 >  if(datasetType_ != "data") {
42 >    if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
43 >    //    muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
44 >    //    electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
45 >  }
46 > #ifdef DISPLACED_SUSY
47 >  if (datasetType_ == "signalMC")
48 >    cTauWeight_ = new CTauWeight (stopCTau_.at (0), stopCTau_.at (1), stops_);
49 > #endif
50 >
51  
52    // Construct Cutflow Objects. These store the results of cut decisions and
53    // handle filling cut flow histograms.
# Line 38 | Line 55 | OSUAnalysis::OSUAnalysis (const edm::Par
55    std::vector<TFileDirectory> directories;
56  
57    //always get vertex collection so we can assign the primary vertex in the event
58 <  
59 < objectsToGet.push_back("primaryvertexs");
60 <  //always make the plot of number of primary verticex (to check pile-up reweighting)
58 >  objectsToGet.push_back("primaryvertexs");
59 >
60 >  //always make the plot of number of primary vertices (to check pile-up reweighting)
61    objectsToPlot.push_back("primaryvertexs");
62  
63    //always get the MC particles to do GEN-matching
# Line 54 | Line 71 | objectsToGet.push_back("primaryvertexs")
71  
72      string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
73      if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
74 <    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
75 <      objectsToGet.push_back(tempInputCollection);
74 >    if(tempInputCollection == "event-track pairs")   tempInputCollection = "track-event pairs";
75 >    if(tempInputCollection == "secondary muon-muon pairs")   tempInputCollection = "muon-secondary muon pairs";
76 >    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
77 >      if(tempInputCollection.find("secondary")!=std::string::npos){//secondary object
78 >        int spaceIndex = tempInputCollection.find(" ");
79 >        int secondWordLength = tempInputCollection.size() - spaceIndex;
80 >        objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
81 >      }
82 >      else{
83 >        objectsToGet.push_back(tempInputCollection);
84 >      }
85        objectsToPlot.push_back(tempInputCollection);
86        objectsToCut.push_back(tempInputCollection);
87      }
88 <    else{//pair of objects, need to add them both to the things to objectsToGet
89 <      int dashIndex = tempInputCollection.find("-");
90 <      int spaceIndex = tempInputCollection.find(" ");
91 <      int secondWordLength = spaceIndex - dashIndex;
92 <      objectsToGet.push_back(tempInputCollection);
67 <      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
68 <      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
69 <      objectsToPlot.push_back(tempInputCollection);
70 <      objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
71 <      objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
88 >    else{//pair of objects, need to add the pair and the individual objects to the lists of things to Get/Plot/Cut
89 >      string obj1;
90 >      string obj2;
91 >      getTwoObjs(tempInputCollection, obj1, obj2);
92 >      string obj2ToGet = getObjToGet(obj2);  
93        objectsToCut.push_back(tempInputCollection);
94 <      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
95 <      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
94 >      objectsToCut.push_back(obj1);
95 >      objectsToCut.push_back(obj2);  
96 >      objectsToPlot.push_back(tempInputCollection);
97 >      objectsToPlot.push_back(obj1);
98 >      objectsToPlot.push_back(obj2);  
99 >      objectsToGet.push_back(tempInputCollection);
100 >      objectsToGet.push_back(obj1);
101 >      objectsToGet.push_back(obj2ToGet);
102  
103 <      }
103 >    }
104  
105      vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
106 <    
106 >
107      for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
108  
109        histogram tempHistogram;
# Line 85 | Line 112 | objectsToGet.push_back("primaryvertexs")
112        tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
113        tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
114        tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
115 <      
115 >
116        histograms.push_back(tempHistogram);
117  
118      }
# Line 100 | Line 127 | objectsToGet.push_back("primaryvertexs")
127    for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
128  
129      string currentObject = objectsToPlot.at(currentObjectIndex);
130 <    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
130 >    if(currentObject != "muons" && currentObject != "secondary muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
131  
132      histogram tempIdHisto;
133      histogram tempMomIdHisto;
134      histogram tempGmaIdHisto;
135 +    histogram tempIdVsMomIdHisto;
136  
137      tempIdHisto.inputCollection = currentObject;
138      tempMomIdHisto.inputCollection = currentObject;
139      tempGmaIdHisto.inputCollection = currentObject;
140 +    tempIdVsMomIdHisto.inputCollection = currentObject;
141 +
142 +    if(currentObject == "secondary muons") currentObject = "secondaryMuons";
143  
144      currentObject = currentObject.substr(0, currentObject.size()-1);
145      tempIdHisto.name = currentObject+"GenMatchId";
146      tempMomIdHisto.name = currentObject+"GenMatchMotherId";
147      tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
148 +    tempIdVsMomIdHisto.name = currentObject+"GenMatchIdVsMotherId";
149  
150      currentObject.at(0) = toupper(currentObject.at(0));
151      tempIdHisto.title = currentObject+" Gen-matched Particle";
152      tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
153      tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
154 +    tempIdVsMomIdHisto.title = currentObject+" Gen-matched Particle's Mother vs. Particle;Particle;Mother";
155  
156      int maxNum = 24;
157      vector<double> binVector;
# Line 132 | Line 165 | objectsToGet.push_back("primaryvertexs")
165      tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
166      tempGmaIdHisto.bins = binVector;
167      tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
168 +    binVector.push_back(maxNum);
169 +    binVector.push_back(0);
170 +    binVector.push_back(maxNum);
171 +    tempIdVsMomIdHisto.bins = binVector;
172 +    tempIdVsMomIdHisto.inputVariables.push_back("genMatchedId");
173 +    tempIdVsMomIdHisto.inputVariables.push_back("genMatchedMotherIdReverse");
174  
175      histograms.push_back(tempIdHisto);
176      histograms.push_back(tempMomIdHisto);
177      histograms.push_back(tempGmaIdHisto);
178 <    
178 >    histograms.push_back(tempIdVsMomIdHisto);
179    }
180  
181  
143
182    channel tempChannel;
183    //loop over all channels (event selections)
184    for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
# Line 177 | Line 215 | objectsToGet.push_back("primaryvertexs")
215      twoDHists_.push_back(twoDhistoMap);
216  
217  
218 +
219      //book all histograms included in the configuration
220      for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
221        histogram currentHistogram = histograms.at(currentHistogramIndex);
# Line 184 | Line 223 | objectsToGet.push_back("primaryvertexs")
223        int numInputVariables = currentHistogram.inputVariables.size();
224  
225        if(numBinsElements != 3 && numBinsElements !=6) {
226 <        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
227 <        exit(0);
226 >        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
227 >        exit(0);
228        }
229        else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
230 <        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
231 <        exit(0);
230 >        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
231 >        exit(0);
232        }
233        else if(numBinsElements == 3){
234 <        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));
234 >        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));
235        }
236        else if(numBinsElements == 6){
237 <        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));
237 >        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));
238  
239        }
240  
241  
242        if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
243  
244 < // bin      particle type
245 < // ---      -------------
246 < //  0        unmatched
247 < //  1        u
248 < //  2        d
249 < //  3        s
250 < //  4        c
251 < //  5        b
252 < //  6        t
253 < //  7        e
254 < //  8        mu
255 < //  9        tau
256 < // 10        nu
257 < // 11        g
258 < // 12        gamma
259 < // 13        Z
260 < // 14        W
261 < // 15        light meson
262 < // 16        K meson
263 < // 17        D meson
264 < // 18        B meson
265 < // 19        light baryon
266 < // 20        strange baryon
267 < // 21        charm baryon
268 < // 22        bottom baryon
269 < // 23        other
244 >      // bin      particle type
245 >      // ---      -------------
246 >      //  0        unmatched
247 >      //  1        u
248 >      //  2        d
249 >      //  3        s
250 >      //  4        c
251 >      //  5        b
252 >      //  6        t
253 >      //  7        e
254 >      //  8        mu
255 >      //  9        tau
256 >      // 10        nu
257 >      // 11        g
258 >      // 12        gamma
259 >      // 13        Z
260 >      // 14        W
261 >      // 15        light meson
262 >      // 16        K meson
263 >      // 17        D meson
264 >      // 18        B meson
265 >      // 19        light baryon
266 >      // 20        strange baryon
267 >      // 21        charm baryon
268 >      // 22        bottom baryon
269 >      // 23        other
270  
271        vector<TString> labelArray;
272        labelArray.push_back("unmatched");
# Line 256 | Line 295 | objectsToGet.push_back("primaryvertexs")
295        labelArray.push_back("other");
296  
297        for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
298 <        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
298 >        if(currentHistogram.name.find("GenMatchIdVsMotherId")==std::string::npos) {
299 >          oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
300 >        }
301 >        else {
302 >          twoDHists_.at(currentChannel)[currentHistogram.name]->GetYaxis()->SetBinLabel(bin+1,labelArray.at(currentHistogram.bins.at(0)-bin-1));
303 >          twoDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
304 >        }
305 >      }
306 >      if(currentHistogram.name.find("GenMatchIdVsMotherId")!=std::string::npos) {
307 >        twoDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->CenterTitle();
308 >        twoDHists_.at(currentChannel)[currentHistogram.name]->GetYaxis()->CenterTitle();
309        }
310 +
311      }
312 <    //book a histogram for the number of each object type to be plotted
312 >
313 >    // Book a histogram for the number of each object type to be plotted.  
314 >    // Name of objectToPlot here must match the name specified in OSUAnalysis::analyze().  
315      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
316        string currentObject = objectsToPlot.at(currentObjectIndex);
317        int maxNum = 10;
318        if(currentObject == "mcparticles") maxNum = 50;
319        else if(currentObject == "primaryvertexs") maxNum = 50;
268      else if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
269      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
270      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
320  
321 <      currentObject.at(0) = toupper(currentObject.at(0));
321 >      if(currentObject == "muon-muon pairs")                currentObject = "dimuonPairs";
322 >      else if(currentObject == "electron-electron pairs")   currentObject = "dielectronPairs";
323 >      else if(currentObject == "electron-muon pairs")       currentObject = "electronMuonPairs";
324 >      else if(currentObject == "track-event pairs")         currentObject = "trackEventPairs";
325 >      else if(currentObject == "electron-track pairs")      currentObject = "electronTrackPairs";        
326 >      else if(currentObject == "muon-track pairs")          currentObject = "muonTrackPairs";    
327 >      else if(currentObject == "muon-tau pairs")            currentObject = "muonTauPairs";      
328 >      else if(currentObject == "tau-tau pairs")             currentObject = "ditauPairs";
329 >      else if(currentObject == "tau-track pairs")           currentObject = "tauTrackPairs";
330 >      else if(currentObject == "muon-secondary muon pairs") currentObject = "muonSecondaryMuonPairs";
331 >      else if(currentObject == "secondary muons")           currentObject = "secondaryMuons";
332 >
333 >      currentObject.at(0) = toupper(currentObject.at(0));  
334        string histoName = "num" + currentObject;
335  
336        if(histoName == "numPrimaryvertexs"){
337 <        string newHistoName = histoName + "BeforePileupCorrection";
338 <        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);    
339 <        newHistoName = histoName + "AfterPileupCorrection";
340 <        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
337 >        string newHistoName = histoName + "BeforePileupCorrection";
338 >        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
339 >        newHistoName = histoName + "AfterPileupCorrection";
340 >        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
341        }
342        else
343 <        oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
343 >        oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
344      }
345  
346  
286    
347  
348  
349      //get list of cuts for this channel
# Line 294 | Line 354 | objectsToGet.push_back("primaryvertexs")
354        cut tempCut;
355        //store input collection for cut
356        string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
357 +      if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
358 +      if(tempInputCollection == "event-track pairs")   tempInputCollection = "track-event pairs";
359 +      if(tempInputCollection == "secondary muon-muon pairs")   tempInputCollection = "muon-secondary muon pairs";
360        tempCut.inputCollection = tempInputCollection;
361 <      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
362 <        objectsToGet.push_back(tempInputCollection);
361 >      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
362 >        if(tempInputCollection.find("secondary")!=std::string::npos){//secondary object
363 >          int spaceIndex = tempInputCollection.find(" ");
364 >          int secondWordLength = tempInputCollection.size() - spaceIndex;
365 >          objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
366 >        }
367 >        else{
368 >          objectsToGet.push_back(tempInputCollection);
369 >        }
370          objectsToCut.push_back(tempInputCollection);
371        }
372 <      else{//pair of objects, need to add them both to the things to objectsToGet
373 <        int dashIndex = tempInputCollection.find("-");
374 <        int spaceIndex = tempInputCollection.find(" ");
375 <        int secondWordLength = spaceIndex - dashIndex;
376 <        objectsToGet.push_back(tempInputCollection);
377 <        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
378 <        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
379 <        objectsToCut.push_back(tempInputCollection);
380 <        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
381 <        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
372 >      else{//pair of objects, need to add them both to objectsToGet
373 >        string obj1;
374 >        string obj2;
375 >        getTwoObjs(tempInputCollection, obj1, obj2);
376 >        string obj2ToGet = getObjToGet(obj2);  
377 >        objectsToCut.push_back(tempInputCollection);
378 >        objectsToCut.push_back(obj1);
379 >        objectsToCut.push_back(obj2);  
380 >        objectsToGet.push_back(tempInputCollection);
381 >        objectsToGet.push_back(obj1);
382 >        objectsToGet.push_back(obj2ToGet);
383  
384        }
385  
# Line 318 | Line 389 | objectsToGet.push_back("primaryvertexs")
389        string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
390        std::vector<string> cutStringVector = splitString(cutString);
391        if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
392 <        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
393 <        exit(0);
392 >        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
393 >        exit(0);
394        }
395        tempCut.numSubcuts = (cutStringVector.size()+1)/4;
396        for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
397 <        int indexOffset = 4 * subcutIndex;
398 <        string currentVariableString = cutStringVector.at(indexOffset);
399 <        if(currentVariableString.find("(")==std::string::npos){
400 <          tempCut.functions.push_back("");//no function was specified
401 <          tempCut.variables.push_back(currentVariableString);// variable to cut on
402 <        }
403 <        else{
404 <          tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
405 <          string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
406 <          tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
407 <        }
408 <        tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
409 <        tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
410 <        if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
397 >        int indexOffset = 4 * subcutIndex;
398 >        string currentVariableString = cutStringVector.at(indexOffset);
399 >        if(currentVariableString.find("(")==std::string::npos){
400 >          tempCut.functions.push_back("");//no function was specified
401 >          tempCut.variables.push_back(currentVariableString);// variable to cut on
402 >        }
403 >        else{
404 >          tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
405 >          string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
406 >          tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
407 >        }
408 >        tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
409 >        tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
410 >        if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
411        }
412  
413        //get number of objects required to pass cut for event to pass
414        string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
415        std::vector<string> numberRequiredVector = splitString(numberRequiredString);
416        if(numberRequiredVector.size()!=2){
417 <        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
417 >        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
418          exit(0);
419        }
420  
# Line 365 | Line 436 | objectsToGet.push_back("primaryvertexs")
436        }
437        tempCut.name = tempCutName;
438  
368
439        tempChannel.cuts.push_back(tempCut);
440  
441  
372
442      }//end loop over cuts
443  
444      channels.push_back(tempChannel);
# Line 445 | Line 514 | OSUAnalysis::analyze (const edm::Event &
514    if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
515      event.getByLabel (superclusters_, superclusters);
516  
517 <
518 <  //get pile-up event weight
519 <  double puScaleFactor = 1.00;
520 <  if(datasetType_ != "data"){
521 <    puScaleFactor = puWeight_->at (events->at (0).numTruePV);
522 <    //    cout << puScaleFactor << endl;
517 >  if (useTrackCaloRhoCorr_) {  
518 >    // Used only for pile-up correction of by-hand calculation of isolation energy.  
519 >    // This rho collection is not available in all BEANs.  
520 >    // For description of rho values for different jet reconstruction algorithms, see
521 >    // https://twiki.cern.ch/twiki/bin/view/CMS/JetAlgorithms#Algorithms  
522 >    event.getByLabel ("kt6CaloJets","rho", rhokt6CaloJetsHandle_);
523    }
524  
525  
526 +  //get pile-up event weight
527 +  double scaleFactor = 1.00;
528 +  if(doPileupReweighting_ && datasetType_ != "data")
529 +    scaleFactor = puWeight_->at (events->at (0).numTruePV);
530 +
531 +  cTauScaleFactor_ = 1.0;
532 + #ifdef DISPLACED_SUSY
533 +  if (datasetType_ == "signalMC")
534 +    cTauScaleFactor_ = cTauWeight_->at (event);
535 + #endif
536 +  scaleFactor *= cTauScaleFactor_;
537  
538    //loop over all channels
539  
# Line 461 | Line 541 | OSUAnalysis::analyze (const edm::Event &
541      channel currentChannel = channels.at(currentChannelIndex);
542  
543      flagMap individualFlags;
464    flagMap cumulativeFlags;
544      counterMap passingCounter;
545 +    cumulativeFlags.clear ();
546  
547      bool triggerDecision = true;
548      if(currentChannel.triggers.size() != 0){  //triggers specified
# Line 480 | Line 560 | OSUAnalysis::analyze (const edm::Event &
560        for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
561          string currentObject = objectsToCut.at(currentObjectIndex);
562  
563 <        //initialize maps to get ready to set cuts
563 >        //initialize maps to get ready to set cuts
564          individualFlags[currentObject].push_back (vector<bool> ());
565          cumulativeFlags[currentObject].push_back (vector<bool> ());
566  
# Line 488 | Line 568 | OSUAnalysis::analyze (const edm::Event &
568        for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
569          string currentObject = objectsToCut.at(currentObjectIndex);
570  
571 <        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
571 >        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
572 >
573  
574          if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
575 +
576          else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
577 +
578 +        else if(currentObject == "secondary muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"secondary muons");
579          else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
580          else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
581          else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
# Line 504 | Line 588 | OSUAnalysis::analyze (const edm::Event &
588          else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
589          else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
590  
507        
508        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
509                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
510                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
511                                                                   "muon-muon pairs");
512        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
513                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
514                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
515                                                                           "electron-electron pairs");
516        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
517                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
518                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
519                                                                       "electron-muon pairs");
520
591  
592 +        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
593 +                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
594 +                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
595 +                                                                   "muon-muon pairs");
596 +
597 +        else if(currentObject == "muon-secondary muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
598 +                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
599 +                                                                   cumulativeFlags.at("secondary muons").at(flagsForPairCutsIndex), \
600 +                                                                   "muon-secondary muon pairs");
601 +
602 +        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
603 +                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
604 +                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
605 +                                                                           "electron-electron pairs");
606 +        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
607 +                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
608 +                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
609 +                                                                       "electron-muon pairs");
610 +        else if(currentObject == "track-event pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),events.product(),
611 +                                                                     cumulativeFlags.at("tracks").at(flagsForPairCutsIndex),
612 +                                                                     cumulativeFlags.at("events").at(flagsForPairCutsIndex),
613 +                                                                     "track-event pairs");
614 +        else if(currentObject == "electron-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),tracks.product(),
615 +                                                                        cumulativeFlags.at("electrons").at(flagsForPairCutsIndex),
616 +                                                                        cumulativeFlags.at("tracks").at(flagsForPairCutsIndex),
617 +                                                                        "electron-track pairs");
618 +        else if(currentObject == "muon-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),tracks.product(),
619 +                                                                    cumulativeFlags.at("muons").at(flagsForPairCutsIndex),
620 +                                                                    cumulativeFlags.at("tracks").at(flagsForPairCutsIndex),
621 +                                                                    "muon-track pairs");
622 +        else if(currentObject == "muon-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),taus.product(),
623 +                                                                  cumulativeFlags.at("muons").at(flagsForPairCutsIndex),
624 +                                                                  cumulativeFlags.at("taus").at(flagsForPairCutsIndex),
625 +                                                                  "muon-tau pairs");
626 +        else if(currentObject == "tau-tau pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),taus.product(),
627 +                                                                 cumulativeFlags.at("taus").at(flagsForPairCutsIndex),
628 +                                                                 cumulativeFlags.at("taus").at(flagsForPairCutsIndex),
629 +                                                                 "tau-tau pairs");
630 +        else if(currentObject == "tau-track pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus .product(),tracks.product(),
631 +                                                                 cumulativeFlags.at("taus").at(flagsForPairCutsIndex),
632 +                                                                 cumulativeFlags.at("tracks").at(flagsForPairCutsIndex),
633 +                                                                 "tau-track pairs");
634 +        
635 +        
636        }
637  
638  
# Line 534 | Line 648 | OSUAnalysis::analyze (const edm::Event &
648      eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
649  
650  
537
651      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
652  
653        //loop over all objects and count how many passed the cumulative selection up to this point
# Line 542 | Line 655 | OSUAnalysis::analyze (const edm::Event &
655        int numberPassing = 0;
656  
657        for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
658 <          if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
658 >        if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
659        }
547
660        bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
661        cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
550
662        eventPassedAllCuts = eventPassedAllCuts && cutDecision;
663  
664      }
665  
666 <    cutFlows_.at(currentChannelIndex)->fillCutFlow(puScaleFactor);
666 >    //     if(datasetType_ != "data") {
667 >    //       scaleFactor *= muonSFWeight_->at (chosenMuon ()->eta);
668 >    //       scaleFactor *= electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
669 >    //     }
670  
671 +    cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
672  
673  
559    if(!eventPassedAllCuts)continue;
560
674  
675 <    
675 >    if(!eventPassedAllCuts)continue;
676  
677 <    //set position of primary vertex in event, in order to calculate quantities relative to it
678 <    if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
679 <      vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
680 <      for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
681 <        if(!vertexFlags.at(vertexIndex)) continue;
682 <        chosenPrimaryVertex = & primaryvertexs->at(vertexIndex);
683 <        break;
571 <      }
572 <    }
573 <    else {
574 <      chosenPrimaryVertex = & primaryvertexs->at(0);
677 >    if (printEventInfo_) {
678 >      // Write information about event to screen, for testing purposes.  
679 >      cout << "Event passed all cuts in channel " <<  currentChannel.name
680 >           << ": run="  << events->at(0).run
681 >           << "  lumi=" << events->at(0).lumi
682 >           << "  event=" << events->at(0).evt
683 >           << endl;  
684      }
685  
686  
578
687      //filling histograms
688      for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
689        histogram currentHistogram = histograms.at(histogramIndex);
690  
691        if(currentHistogram.inputVariables.size() == 1){
692 <        TH1D* histo;
693 <        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
586 <
587 <
692 >        TH1D* histo;
693 >        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
694  
695 <        if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
696 <        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
697 <        else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),\
698 <                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(),\
699 <                                                                                       cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
700 <        else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),puScaleFactor);
701 <        else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
702 <                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
703 <                                                                                               cumulativeFlags.at("electron-electron pairs").back(),puScaleFactor);
704 <        else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
705 <                                                                                              cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
706 <                                                                                              cumulativeFlags.at("electron-muon pairs").back(),puScaleFactor);
707 <        else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),puScaleFactor);
708 <        else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),puScaleFactor);
709 <        else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),puScaleFactor);
710 <        else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),puScaleFactor);
711 <        else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),puScaleFactor);
712 <        else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),puScaleFactor);
713 <        else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),puScaleFactor);
714 <        else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),puScaleFactor);
715 <        else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),puScaleFactor);
716 <        else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),puScaleFactor);
695 >        if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
696 >        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
697 >        else if(currentHistogram.inputCollection == "secondary muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").back(),scaleFactor);
698 >        else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
699 >                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
700 >                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
701 >        else if(currentHistogram.inputCollection == "muon-secondary muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
702 >                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("secondary muons").back(), \
703 >                                                                                       cumulativeFlags.at("muon-secondary muon pairs").back(),scaleFactor);
704 >        else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
705 >        else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
706 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
707 >                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
708 >        else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
709 >                                                                                           cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
710 >                                                                                           cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
711 >        else if(currentHistogram.inputCollection == "electron-track pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),tracks.product(),        
712 >                                                                                            cumulativeFlags.at("electrons").back(),cumulativeFlags.at("tracks").back(),  
713 >                                                                                            cumulativeFlags.at("electron-track pairs").back(),scaleFactor);      
714 >        else if(currentHistogram.inputCollection == "muon-track pairs") fill1DHistogram(histo,currentHistogram, muons.product(),tracks.product(),        
715 >                                                                                        cumulativeFlags.at("muons").back(),cumulativeFlags.at("tracks").back(),  
716 >                                                                                        cumulativeFlags.at("muon-track pairs").back(),scaleFactor);      
717 >        else if(currentHistogram.inputCollection == "muon-tau pairs") fill1DHistogram(histo,currentHistogram, muons.product(),taus.product(),    
718 >                                                                                      cumulativeFlags.at("muons").back(),cumulativeFlags.at("taus").back(),      
719 >                                                                                      cumulativeFlags.at("muon-tau pairs").back(),scaleFactor);  
720 >        else if(currentHistogram.inputCollection == "tau-tau pairs") fill1DHistogram(histo,currentHistogram, taus.product(),taus.product(),      
721 >                                                                                     cumulativeFlags.at("taus").back(),cumulativeFlags.at("taus").back(),        
722 >                                                                                     cumulativeFlags.at("tau-tau pairs").back(),scaleFactor);    
723 >        else if(currentHistogram.inputCollection == "tau-track pairs") fill1DHistogram(histo,currentHistogram, taus.product(),tracks.product(),
724 >                                                                                     cumulativeFlags.at("taus").back(),cumulativeFlags.at("tracks").back(),
725 >                                                                                     cumulativeFlags.at("tau-track pairs").back(),scaleFactor);
726 >
727 >        else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
728 >        else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
729 >        else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
730 >        else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
731 >        else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
732 >        else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
733 >        else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
734 >        else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
735 >        else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
736 >        else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
737        }
738        else if(currentHistogram.inputVariables.size() == 2){
739 <        TH2D* histo;
740 <        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
739 >        TH2D* histo;
740 >        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
741  
742 +        if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
743 +        else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
744 +        else if(currentHistogram.inputCollection == "secondary muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").back(),scaleFactor);
745 +        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
746 +                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
747 +                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
748 +        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
749 +                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("secondary muons").back(), \
750 +                                                                                       cumulativeFlags.at("muon-secondary muon pairs").back(),scaleFactor);
751 +        else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
752 +        else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
753 +                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
754 +                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
755 +        else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
756 +                                                                                           cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
757 +                                                                                           cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
758 +        else if(currentHistogram.inputCollection == "electron-track pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),tracks.product(),        
759 +                                                                                            cumulativeFlags.at("electrons").back(),cumulativeFlags.at("tracks").back(),          
760 +                                                                                            cumulativeFlags.at("electron-track pairs").back(),scaleFactor);      
761 +        else if(currentHistogram.inputCollection == "muon-track pairs") fill2DHistogram(histo,currentHistogram,muons.product(),tracks.product(),        
762 +                                                                                        cumulativeFlags.at("muons").back(),cumulativeFlags.at("tracks").back(),          
763 +                                                                                        cumulativeFlags.at("muon-track pairs").back(),scaleFactor);      
764 +        else if(currentHistogram.inputCollection == "muon-tau pairs") fill2DHistogram(histo,currentHistogram,muons.product(),taus.product(),    
765 +                                                                                      cumulativeFlags.at("muons").back(),cumulativeFlags.at("taus").back(),      
766 +                                                                                      cumulativeFlags.at("muon-tau pairs").back(),scaleFactor);  
767 +        else if(currentHistogram.inputCollection == "tau-tau pairs") fill2DHistogram(histo,currentHistogram,taus.product(),taus.product(),      
768 +                                                                                     cumulativeFlags.at("taus").back(),cumulativeFlags.at("taus").back(),        
769 +                                                                                     cumulativeFlags.at("tau-tau pairs").back(),scaleFactor);    
770 +        else if(currentHistogram.inputCollection == "tau-track pairs") fill2DHistogram(histo,currentHistogram,taus.product(),tracks.product(),
771 +                                                                                     cumulativeFlags.at("taus").back(),cumulativeFlags.at("tracks").back(),
772 +                                                                                     cumulativeFlags.at("tau-track pairs").back(),scaleFactor);
773 +        else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
774 +        else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
775 +        else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
776 +        else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
777 +        else if(currentHistogram.inputCollection == "track-event pairs") fill2DHistogram(histo,currentHistogram,tracks.product(),events.product(),
778 +                                                                                         cumulativeFlags.at("tracks").back(),cumulativeFlags.at("events").back(),
779 +                                                                                         cumulativeFlags.at("track-event pairs").back(),scaleFactor);  
780 +        else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
781 +        else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
782 +        else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
783 +        else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
784 +        else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
785 +        else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
786 +      }
787 +    }
788  
789  
618        if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
619        else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
620        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
621                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
622                                                                                       cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
623        else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),puScaleFactor);
624        else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
625                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
626                                                                                               cumulativeFlags.at("electron-electron pairs").back(),puScaleFactor);
627        else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
628                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
629                                                                                               cumulativeFlags.at("electron-muon pairs").back(),puScaleFactor);
630        else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),puScaleFactor);
631        else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),puScaleFactor);
632        else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),puScaleFactor);
633        else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),puScaleFactor);
634        else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),puScaleFactor);
635        else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),puScaleFactor);
636        else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),puScaleFactor);
637        else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),puScaleFactor);
638        else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),puScaleFactor);
639        else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),puScaleFactor);
640      }
641    }
790  
791      //fills histograms with the sizes of collections
792      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
793 +
794        string currentObject = objectsToPlot.at(currentObjectIndex);
795  
796 <      string objectToPlot = "";
796 >      string objectToPlot = "";  
797  
798 +      // Name of objectToPlot here must match the name specified in OSUAnalysis::OSUAnalysis().  
799        if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
800        else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
801 <      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
801 >      else if(currentObject == "electron-muon pairs")     objectToPlot = "electronMuonPairs";
802 >      else if(currentObject == "electron-track pairs")    objectToPlot = "electronTrackPairs";
803 >      else if(currentObject == "muon-track pairs")        objectToPlot = "muonTrackPairs";      
804 >      else if(currentObject == "muon-tau pairs")          objectToPlot = "muonTauPairs";        
805 >      else if(currentObject == "tau-tau pairs")           objectToPlot = "ditauPairs";
806 >      else if(currentObject == "tau-track pairs")         objectToPlot = "tauTrackPairs";
807 >      else if(currentObject == "track-event pairs")       objectToPlot = "trackEventPairs";  
808 >      else if(currentObject == "muon-secondary muon pairs")       objectToPlot = "muonSecondaryMuonPairs";  
809 >      else if(currentObject == "secondary muons")         objectToPlot = "secondaryMuons";  
810        else objectToPlot = currentObject;
811 +
812        string tempCurrentObject = objectToPlot;
813 <      tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
813 >      tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));  
814        string histoName = "num" + tempCurrentObject;
815  
657
658
659
816        //set position of primary vertex in event, in order to calculate quantities relative to it
817        if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
818 <        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
818 >        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
819          int numToPlot = 0;
820 <        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
821 <          if(lastCutFlags.at(currentFlag)) numToPlot++;
822 <        }
823 <        if(objectToPlot == "primaryvertexs"){
824 <          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
825 <          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
820 >        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
821 >          if(lastCutFlags.at(currentFlag)) numToPlot++;
822 >        }
823 >        if(objectToPlot == "primaryvertexs"){
824 >          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
825 >          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
826 >        }
827 >        else {
828 >          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,scaleFactor);
829          }
671        else
672          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,puScaleFactor);        
673      }
674      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
675      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
676      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
677      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
678      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
679      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
680      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
681      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
682      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
683      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
684      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
685      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
686      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor);
687      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor);
688      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor);
689      else if(objectToPlot == "primaryvertexs"){
690        oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
691        oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
830        }
831 +      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),scaleFactor);
832 +      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),scaleFactor);
833 +      else if(objectToPlot == "secondaryMuons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),scaleFactor);
834 +      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,scaleFactor);
835 +      else if(objectToPlot == "muonSecondaryMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,scaleFactor);
836 +      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),scaleFactor);
837 +      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,scaleFactor);
838 +      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),scaleFactor);
839 +      else if(objectToPlot == "electronTrackPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*tracks->size(),scaleFactor);
840 +      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),scaleFactor);
841 +      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),scaleFactor);
842 +      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),scaleFactor);
843 +      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),scaleFactor);
844 +      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),scaleFactor);
845 +      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),scaleFactor);
846 +      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),scaleFactor);
847 +      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),scaleFactor);
848 +      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),scaleFactor);
849 +      else if(objectToPlot == "primaryvertexs"){
850 +        oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
851 +        oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
852 +      }
853  
854 <    }
695 <
696 <
697 <
854 >    } // end for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
855  
856    } //end loop over channel
857  
858 <  masterCutFlow_->fillCutFlow(puScaleFactor);
858 >  masterCutFlow_->fillCutFlow(scaleFactor);
859  
860 + } // end void OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
861  
862  
705 }
706
863  
864   bool
865   OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
# Line 750 | Line 906 | OSUAnalysis::splitString (string inputSt
906  
907   }
908  
909 +
910 + void OSUAnalysis::getTwoObjs(string tempInputCollection, string& obj1, string& obj2) {
911 +  // Set two object strings from the tempInputCollection string,
912 +  // For example, if tempInputCollection is "electron-muon pairs",
913 +  // then obj1 = "electrons" and obj2 = "muons".  
914 +  // Note that the objects have an "s" appended.  
915 +
916 +  int dashIndex = tempInputCollection.find("-");
917 +  int spaceIndex = tempInputCollection.find_last_of(" ");
918 +  int secondWordLength = spaceIndex - dashIndex;
919 +  obj1 = tempInputCollection.substr(0,dashIndex) + "s";  
920 +  obj2 = tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s";
921 +  
922 + }
923 +
924 +
925 + string OSUAnalysis::getObjToGet(string obj) {
926 +  // Return the string corresponding to the object to get for the given obj string.
927 +  // Right now this only handles the case in which obj contains "secondary",
928 +  // e.g, "secondary muons".  
929 +  // Note that "s" is NOT appended.  
930 +  
931 +  if (obj.find("secondary")==std::string::npos) return obj;  // "secondary" is not found  
932 +  int firstSpaceIndex = obj.find_first_of(" ");  
933 +  return obj.substr(firstSpaceIndex+1,obj.length()-1);
934 +
935 + }  
936 +
937 +
938 +
939   double
940   OSUAnalysis::valueLookup (const BNjet* object, string variable, string function){
941  
# Line 1036 | Line 1222 | OSUAnalysis::valueLookup (const BNmuon*
1222    else if(variable == "time_ndof") value = object->time_ndof;
1223  
1224    //user-defined variables
1225 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1226 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1225 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1226 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1227    else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1228 <  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1228 >  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1229    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1230 +  else if(variable == "metMT") {
1231 +    const BNmet *met = chosenMET ();
1232 +    if (met)
1233 +      {
1234 +        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1235 +          p2 (met->px, met->py, 0.0, met->pt);
1236 +
1237 +        value = (p1 + p2).Mt ();
1238 +      }
1239 +    else
1240 +      value = -999;
1241 +  }
1242  
1243  
1244  
# Line 1151 | Line 1349 | OSUAnalysis::valueLookup (const BNmuon*
1349  
1350  
1351    else if(variable == "tightID") {
1352 <    value = object->isGlobalMuon > 0            \
1353 <      && object->isPFMuon > 0                   \
1354 <      && object->normalizedChi2 < 10            \
1355 <      && object->numberOfValidMuonHits > 0      \
1356 <      && object->numberOfMatchedStations > 1    \
1357 <      && fabs(object->correctedD0Vertex) < 0.2  \
1358 <      && fabs(object->correctedDZ) < 0.5        \
1359 <      && object->numberOfValidPixelHits > 0             \
1352 >    value = object->isGlobalMuon > 0                \
1353 >      && object->isPFMuon > 0                        \
1354 >      && object->normalizedChi2 < 10                \
1355 >                                  && object->numberOfValidMuonHits > 0        \
1356 >      && object->numberOfMatchedStations > 1        \
1357 >      && fabs(object->correctedD0Vertex) < 0.2        \
1358 >      && fabs(object->correctedDZ) < 0.5        \
1359 >      && object->numberOfValidPixelHits > 0                \
1360        && object->numberOfLayersWithMeasurement > 5;
1361    }
1362    else if(variable == "tightIDdisplaced"){
1363 <    value = object->isGlobalMuon > 0            \
1364 <      && object->normalizedChi2 < 10            \
1365 <      && object->numberOfValidMuonHits > 0      \
1366 <      && object->numberOfMatchedStations > 1    \
1367 <      && object->numberOfValidPixelHits > 0     \
1363 >    value = object->isGlobalMuon > 0                \
1364 >      && object->normalizedChi2 < 10                \
1365 >                                  && object->numberOfValidMuonHits > 0        \
1366 >      && object->numberOfMatchedStations > 1        \
1367 >      && object->numberOfValidPixelHits > 0        \
1368        && object->numberOfLayersWithMeasurement > 5;
1369    }
1370  
1371 +  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1372  
1373 +  else if(variable == "genMatchedPdgId"){
1374 +    int index = getGenMatchedParticleIndex(object);
1375 +    if(index == -1) value = 0;
1376 +    else value = mcparticles->at(index).id;
1377 +  }
1378  
1379    else if(variable == "genMatchedId"){
1380      int index = getGenMatchedParticleIndex(object);
# Line 1182 | Line 1386 | OSUAnalysis::valueLookup (const BNmuon*
1386      if(index == -1) value = 0;
1387      else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1388    }
1389 +  else if(variable == "genMatchedMotherIdReverse"){
1390 +    int index = getGenMatchedParticleIndex(object);
1391 +    if(index == -1) value = 23;
1392 +    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1393 +  }
1394    else if(variable == "genMatchedGrandmotherId"){
1395      int index = getGenMatchedParticleIndex(object);
1396      if(index == -1) value = 0;
# Line 1365 | Line 1574 | OSUAnalysis::valueLookup (const BNelectr
1574    else if(variable == "passConvVeto") value = object->passConvVeto;
1575  
1576    //user-defined variables
1577 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1578 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1577 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1578 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1579    else if(variable == "detIso") value = (object->trackIso) / object->pt;
1580    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1581 +  else if(variable == "metMT") {
1582 +    const BNmet *met = chosenMET ();
1583 +    if (met)
1584 +      {
1585 +        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1586 +          p2 (met->px, met->py, 0.0, met->pt);
1587  
1588 +        value = (p1 + p2).Mt ();
1589 +      }
1590 +    else
1591 +      value = -999;
1592 +  }
1593  
1594    else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
1595      if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
# Line 1380 | Line 1600 | OSUAnalysis::valueLookup (const BNelectr
1600      else value = -999;
1601    }
1602  
1383
1603    else if(variable == "correctedD0VertexInEBPlus"){
1604      if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1605      else value = -999;
# Line 1433 | Line 1652 | OSUAnalysis::valueLookup (const BNelectr
1652      else value = -999;
1653    }
1654  
1655 +  else if(variable == "tightID"){
1656 +    if (object->isEB)
1657 +      {
1658 +        value = fabs(object->delEtaIn) < 0.004 \
1659 +          && fabs (object->delPhiIn) < 0.03 \
1660 +          && object->scSigmaIEtaIEta < 0.01 \
1661 +          && object->hadOverEm < 0.12 \
1662 +          && fabs (object->correctedD0Vertex) < 0.02 \
1663 +          && fabs (object->correctedDZ) < 0.1 \
1664 +          && object->absInvEMinusInvPin < 0.05 \
1665 +          && object->passConvVeto;
1666 +      }
1667 +    else
1668 +      {
1669 +        value = fabs(object->delEtaIn) < 0.005 \
1670 +          && fabs (object->delPhiIn) < 0.02 \
1671 +          && object->scSigmaIEtaIEta < 0.03 \
1672 +          && object->hadOverEm < 0.10 \
1673 +          && fabs (object->correctedD0Vertex) < 0.02 \
1674 +          && fabs (object->correctedDZ) < 0.1 \
1675 +          && object->absInvEMinusInvPin < 0.05 \
1676 +          && object->passConvVeto;
1677 +      }
1678 +  }
1679  
1680    else if(variable == "correctedD0VertexInEBPositiveCharge"){
1681      if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
# Line 1490 | Line 1733 | OSUAnalysis::valueLookup (const BNelectr
1733    else if(variable == "tightIDdisplaced"){
1734      if (object->isEB)
1735        {
1736 <        value = fabs(object->delEtaIn) < 0.004 \
1737 <          && fabs (object->delPhiIn) < 0.03 \
1736 >        value = fabs(object->delEtaIn) < 0.004 \
1737 >          && fabs (object->delPhiIn) < 0.03 \
1738            && object->scSigmaIEtaIEta < 0.01 \
1739            && object->hadOverEm < 0.12 \
1740 <          && object->absInvEMinusInvPin < 0.05;
1740 >          && object->absInvEMinusInvPin < 0.05;
1741        }
1742      else
1743        {
1744 <        value = fabs (object->delEtaIn) < 0.005 \
1745 <          && fabs (object->delPhiIn) < 0.02 \
1744 >        value = fabs (object->delEtaIn) < 0.005 \
1745 >          && fabs (object->delPhiIn) < 0.02 \
1746            && object->scSigmaIEtaIEta < 0.03 \
1747            && object->hadOverEm < 0.10 \
1748 <          && object->absInvEMinusInvPin < 0.05;
1748 >          && object->absInvEMinusInvPin < 0.05;
1749        }
1750    }
1751  
1752 +
1753 +  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1754 +
1755 +  else if(variable == "genMatchedPdgId"){
1756 +    int index = getGenMatchedParticleIndex(object);
1757 +    if(index == -1) value = 0;
1758 +    else value = mcparticles->at(index).id;
1759 +  }
1760 +
1761 +
1762    else if(variable == "genMatchedId"){
1763      int index = getGenMatchedParticleIndex(object);
1764      if(index == -1) value = 0;
# Line 1516 | Line 1769 | OSUAnalysis::valueLookup (const BNelectr
1769      if(index == -1) value = 0;
1770      else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1771    }
1772 +  else if(variable == "genMatchedMotherIdReverse"){
1773 +    int index = getGenMatchedParticleIndex(object);
1774 +    if(index == -1) value = 23;
1775 +    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1776 +  }
1777    else if(variable == "genMatchedGrandmotherId"){
1778      int index = getGenMatchedParticleIndex(object);
1779      if(index == -1) value = 0;
# Line 1606 | Line 1864 | OSUAnalysis::valueLookup (const BNevent*
1864    else if(variable == "id1") value = object->id1;
1865    else if(variable == "id2") value = object->id2;
1866    else if(variable == "evt") value = object->evt;
1867 +  else if(variable == "puScaleFactor"){
1868 +    if(datasetType_ != "data")
1869 +      value = puWeight_->at (events->at (0).numTruePV);
1870 +    else
1871 +      value = 1.0;
1872 +  }
1873 +  else if(variable == "muonScaleFactor"){
1874 +    if(datasetType_ != "data")
1875 +      //      value = muonSFWeight_->at (chosenMuon ()->eta);
1876 +      value = 1.0;
1877 +    else
1878 +      value = 1.0;
1879 +  }
1880 +  else if(variable == "electronScaleFactor"){
1881 +    if(datasetType_ != "data")
1882 +      //      value = electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
1883 +      value = 1.0;
1884 +    else
1885 +      value = 1.0;
1886 +  }
1887 +  else if(variable == "cTauScaleFactor")
1888 +    value = cTauScaleFactor_;
1889  
1890    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1891  
# Line 1660 | Line 1940 | OSUAnalysis::valueLookup (const BNtau* o
1940    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1941  
1942  
1943 +
1944 +  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1945 +
1946 +  else if(variable == "genMatchedPdgId"){
1947 +    int index = getGenMatchedParticleIndex(object);
1948 +    if(index == -1) value = 0;
1949 +    else value = mcparticles->at(index).id;
1950 +  }
1951 +
1952    else if(variable == "genMatchedId"){
1953      int index = getGenMatchedParticleIndex(object);
1954      if(index == -1) value = 0;
# Line 1670 | Line 1959 | OSUAnalysis::valueLookup (const BNtau* o
1959      if(index == -1) value = 0;
1960      else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1961    }
1962 +  else if(variable == "genMatchedMotherIdReverse"){
1963 +    int index = getGenMatchedParticleIndex(object);
1964 +    if(index == -1) value = 23;
1965 +    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1966 +  }
1967    else if(variable == "genMatchedGrandmotherId"){
1968      int index = getGenMatchedParticleIndex(object);
1969      if(index == -1) value = 0;
# Line 1766 | Line 2060 | OSUAnalysis::valueLookup (const BNtrack*
2060  
2061    double value = 0.0;
2062    double pMag = sqrt(object->pt * object->pt +
2063 <                         object->pz * object->pz);
2063 >                     object->pz * object->pz);
2064    
2065    if(variable == "pt") value = object->pt;
2066    else if(variable == "px") value = object->px;
# Line 1785 | Line 2079 | OSUAnalysis::valueLookup (const BNtrack*
2079    else if(variable == "charge") value = object->charge;
2080    else if(variable == "numValidHits") value = object->numValidHits;
2081    else if(variable == "isHighPurity") value = object->isHighPurity;
2082 <  //additional BNs info for disappTrks  
2082 >
2083 >
2084 >  //additional BNs info for disappTrks
2085    else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
2086    else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
2087    else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
# Line 1796 | Line 2092 | OSUAnalysis::valueLookup (const BNtrack*
2092    else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
2093    else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
2094    else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
2095 +  
2096 +
2097    //user defined variables
2098    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;
2099    else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
2100 <  else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
2101 <  else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
2102 <  else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
2103 <  else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
2100 >  else if(variable == "caloTotDeltaRp5")            value =  (object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
2101 >  else if(variable == "caloTotDeltaRp5ByP")         value = ((object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
2102 >  else if(variable == "caloTotDeltaRp5RhoCorr")     value = getTrkCaloTotRhoCorr(object);  
2103 >  else if(variable == "caloTotDeltaRp5ByPRhoCorr")  value = getTrkCaloTotRhoCorr(object) / pMag;  
2104 >  else if(variable == "isIso")                      value = getTrkIsIso(object, tracks.product());
2105 >  else if(variable == "isMatchedDeadEcal")          value = getTrkIsMatchedDeadEcal(object);
2106 >  else if(variable == "ptErrorByPt")                value = (object->ptError/object->pt);
2107 >  else if(variable == "ptError")                    value = object->ptError;
2108 >  else if(variable == "ptRes")                      value = getTrkPtRes(object);
2109 >  else if (variable == "d0wrtPV"){      
2110 >    double vx = object->vx - chosenVertex ()->x,        
2111 >      vy = object->vy - chosenVertex ()->y,      
2112 >      px = object->px,  
2113 >      py = object->py,  
2114 >      pt = object->pt;  
2115 >    value = (-vx * py + vy * px) / pt;  
2116 >  }      
2117 >  else if (variable == "dZwrtPV"){      
2118 >    double vx = object->vx - chosenVertex ()->x,        
2119 >      vy = object->vy - chosenVertex ()->y,      
2120 >      vz = object->vz - chosenVertex ()->z,      
2121 >      px = object->px,  
2122 >      py = object->py,  
2123 >      pz = object->pz,  
2124 >      pt = object->pt;  
2125 >    value = vz - (vx * px + vy * py)/pt * (pz/pt);      
2126 >  }    
2127 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2128 >
2129 >  else if(variable == "genMatchedPdgId"){
2130 >    int index = getGenMatchedParticleIndex(object);
2131 >    if(index == -1) value = 0;
2132 >    else value = mcparticles->at(index).id;
2133 >  }
2134  
2135  
2136    else if(variable == "genMatchedId"){
# Line 1815 | Line 2143 | OSUAnalysis::valueLookup (const BNtrack*
2143      if(index == -1) value = 0;
2144      else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2145    }
2146 +  else if(variable == "genMatchedMotherIdReverse"){
2147 +    int index = getGenMatchedParticleIndex(object);
2148 +    if(index == -1) value = 23;
2149 +    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2150 +  }
2151    else if(variable == "genMatchedGrandmotherId"){
2152      int index = getGenMatchedParticleIndex(object);
2153      if(index == -1) value = 0;
# Line 1955 | Line 2288 | OSUAnalysis::valueLookup (const BNmcpart
2288  
2289    //user-defined variables
2290    else if (variable == "d0"){
2291 <    double vx = object->vx - chosenPrimaryVertex->x,
2292 <      vy = object->vy - chosenPrimaryVertex->y,
2291 >    double vx = object->vx - chosenVertex ()->x,
2292 >      vy = object->vy - chosenVertex ()->y,
2293        px = object->px,
2294        py = object->py,
2295        pt = object->pt;
2296      value = (-vx * py + vy * px) / pt;
2297    }
2298    else if (variable == "dz"){
2299 <    double vx = object->vx - chosenPrimaryVertex->x,
2300 <      vy = object->vy - chosenPrimaryVertex->y,
2301 <      vz = object->vz - chosenPrimaryVertex->z,
2299 >    double vx = object->vx - chosenVertex ()->x,
2300 >      vy = object->vy - chosenVertex ()->y,
2301 >      vz = object->vz - chosenVertex ()->z,
2302        px = object->px,
2303        py = object->py,
2304        pz = object->pz,
# Line 1976 | Line 2309 | OSUAnalysis::valueLookup (const BNmcpart
2309      value = sqrt(object->vx*object->vx + object->vy*object->vy);
2310    }
2311    else if(variable == "deltaV0"){
2312 <    value = sqrt((object->vx-chosenPrimaryVertex->x)*(object->vx-chosenPrimaryVertex->x) + (object->vy-chosenPrimaryVertex->y)*(object->vy-chosenPrimaryVertex->y));
2312 >    value = sqrt((object->vx-chosenVertex ()->x)*(object->vx-chosenVertex ()->x) + (object->vy-chosenVertex ()->y)*(object->vy-chosenVertex ()->y));
2313    }
2314    else if (variable == "deltaVx"){
2315 <    value = object->vx - chosenPrimaryVertex->x;
2315 >    value = object->vx - chosenVertex ()->x;
2316    }
2317    else if (variable == "deltaVy"){
2318 <    value = object->vy - chosenPrimaryVertex->y;
2318 >    value = object->vy - chosenVertex ()->y;
2319    }
2320    else if (variable == "deltaVz"){
2321 <    value = object->vz - chosenPrimaryVertex->z;
2321 >    value = object->vz - chosenVertex ()->z;
2322    }
2323  
2324  
# Line 2116 | Line 2449 | OSUAnalysis::valueLookup (const BNphoton
2449    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
2450  
2451  
2452 +
2453 +
2454 +  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2455 +
2456 +  else if(variable == "genMatchedPdgId"){
2457 +    int index = getGenMatchedParticleIndex(object);
2458 +    if(index == -1) value = 0;
2459 +    else value = mcparticles->at(index).id;
2460 +  }
2461 +
2462 +
2463 +
2464    else if(variable == "genMatchedId"){
2465      int index = getGenMatchedParticleIndex(object);
2466      if(index == -1) value = 0;
# Line 2126 | Line 2471 | OSUAnalysis::valueLookup (const BNphoton
2471      if(index == -1) value = 0;
2472      else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2473    }
2474 +  else if(variable == "genMatchedMotherIdReverse"){
2475 +    int index = getGenMatchedParticleIndex(object);
2476 +    if(index == -1) value = 23;
2477 +    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2478 +  }
2479    else if(variable == "genMatchedGrandmotherId"){
2480      int index = getGenMatchedParticleIndex(object);
2481      if(index == -1) value = 0;
# Line 2174 | Line 2524 | OSUAnalysis::valueLookup (const BNmuon*
2524    double value = 0.0;
2525  
2526    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2527 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2528    else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2529    else if(variable == "invMass"){
2530      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2531      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2532 <
2533 <    value = (fourVector1 + fourVector2).M();}
2532 >    value = (fourVector1 + fourVector2).M();
2533 >  }
2534 >  else if(variable == "pt"){
2535 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2536 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2537 >    value = (fourVector1 + fourVector2).Pt();
2538 >  }
2539    else if(variable == "threeDAngle")
2540      {
2541        TVector3 threeVector1(object1->px, object1->py, object1->pz);
2542        TVector3 threeVector2(object2->px, object2->py, object2->pz);
2543        value = (threeVector1.Angle(threeVector2));
2544      }
2545 <
2546 <
2547 <
2545 >  else if(variable == "alpha")
2546 >    {
2547 >      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
2548 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2549 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2550 >      value = (pi-threeVector1.Angle(threeVector2));
2551 >    }
2552    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2553    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2554    else if(variable == "d0Sign"){
# Line 2197 | Line 2557 | OSUAnalysis::valueLookup (const BNmuon*
2557      else if (d0Sign > 0) value = 0.5;
2558      else value = -999;
2559    }
2560 +  else if(variable == "chargeProduct"){
2561 +    value = object1->charge*object2->charge;
2562 +  }
2563    else if(variable == "muon1CorrectedD0Vertex"){
2564      value = object1->correctedD0Vertex;
2565    }
2566    else if(variable == "muon2CorrectedD0Vertex"){
2567      value = object2->correctedD0Vertex;
2568    }
2569 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2570 <
2571 <
2572 <
2573 <
2574 <
2575 <
2569 >  else if(variable == "muon1timeAtIpInOut"){
2570 >    value = object1->timeAtIpInOut;
2571 >  }
2572 >  else if(variable == "muon2timeAtIpInOut"){
2573 >    value = object2->timeAtIpInOut;
2574 >  }
2575 >  else if(variable == "muon1correctedD0")
2576 >    {
2577 >      value = object1->correctedD0;
2578 >    }
2579 >  else if(variable == "muon2correctedD0")
2580 >    {
2581 >      value = object2->correctedD0;
2582 >    }
2583  
2584 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2585  
2586    value = applyFunction(function, value);
2587  
# Line 2223 | Line 2594 | OSUAnalysis::valueLookup (const BNelectr
2594    double value = 0.0;
2595  
2596    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2597 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2598    else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2599    else if(variable == "invMass"){
2600      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2601      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2602 <
2603 <    value = (fourVector1 + fourVector2).M();}
2604 < else if(variable == "threeDAngle")
2602 >    value = (fourVector1 + fourVector2).M();
2603 >  }
2604 >  else if(variable == "pt"){
2605 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2606 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2607 >    value = (fourVector1 + fourVector2).Pt();
2608 >  }
2609 >  else if(variable == "threeDAngle")
2610      {
2611        TVector3 threeVector1(object1->px, object1->py, object1->pz);
2612        TVector3 threeVector2(object2->px, object2->py, object2->pz);
2613        value = (threeVector1.Angle(threeVector2));
2614      }
2238  
2239
2240
2241
2615    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2616    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2244
2245  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2246
2247
2248
2617    else if(variable == "d0Sign"){
2618      double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2619      if(d0Sign < 0) value = -0.5;
2620      else if (d0Sign > 0) value = 0.5;
2621      else value = -999;
2622    }
2623 +  else if(variable == "chargeProduct"){
2624 +    value = object1->charge*object2->charge;
2625 +  }
2626    else if(variable == "electron1CorrectedD0Vertex"){
2627      value = object1->correctedD0Vertex;
2628    }
2629    else if(variable == "electron2CorrectedD0Vertex"){
2630      value = object2->correctedD0Vertex;
2631    }
2632 +  else if(variable == "electron1CorrectedD0"){
2633 +    value = object1->correctedD0;
2634 +  }
2635 +  else if(variable == "electron2CorrectedD0"){
2636 +    value = object2->correctedD0;
2637 +  }
2638  
2639    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2640  
# Line 2272 | Line 2649 | OSUAnalysis::valueLookup (const BNelectr
2649    double value = 0.0;
2650  
2651    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2652 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2653    else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2654    else if(variable == "invMass"){
2655      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2656      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2657 <
2658 <    value = (fourVector1 + fourVector2).M();}
2659 < else if(variable == "threeDAngle")
2657 >    value = (fourVector1 + fourVector2).M();
2658 >  }
2659 >  else if(variable == "pt"){
2660 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2661 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2662 >    value = (fourVector1 + fourVector2).Pt();
2663 >  }
2664 >  else if(variable == "threeDAngle")
2665      {
2666        TVector3 threeVector1(object1->px, object1->py, object1->pz);
2667        TVector3 threeVector2(object2->px, object2->py, object2->pz);
2668        value = (threeVector1.Angle(threeVector2));
2669      }
2287
2288
2670    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2671    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2672    else if(variable == "d0Sign"){
# Line 2294 | Line 2675 | OSUAnalysis::valueLookup (const BNelectr
2675      else if (d0Sign > 0) value = 0.5;
2676      else value = -999;
2677    }
2678 +  else if(variable == "chargeProduct"){
2679 +    value = object1->charge*object2->charge;
2680 +  }
2681    else if(variable == "electronCorrectedD0Vertex"){
2682      value = object1->correctedD0Vertex;
2683    }
# Line 2306 | Line 2690 | OSUAnalysis::valueLookup (const BNelectr
2690    else if(variable == "muonCorrectedD0"){
2691      value = object2->correctedD0;
2692    }
2693 +  else if(variable == "electronDetIso"){
2694 +    value = (object1->trackIso) / object1->pt;
2695 +  }
2696 +  else if(variable == "muonDetIso"){
2697 +    value = (object2->trackIsoDR03) / object2->pt;
2698 +  }
2699 +  else if(variable == "electronRelPFrhoIso"){
2700 +    value = ( object1->chargedHadronIsoDR03 + max(0.0, object1->neutralHadronIsoDR03 + object1->photonIsoDR03 - object1->AEffDr03*object1->rhoPrime) ) / object1->pt;
2701 +  }
2702 +  else if(variable == "muonRelPFdBetaIso"){
2703 +    value = (object2->pfIsoR04SumChargedHadronPt + max(0.0, object2->pfIsoR04SumNeutralHadronEt + object2->pfIsoR04SumPhotonEt - 0.5*object2->pfIsoR04SumPUPt)) / object2->pt;
2704 +  }
2705 +
2706 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2707 +
2708 +  value = applyFunction(function, value);
2709 +
2710 +  return value;
2711 + }
2712 +
2713 +
2714 + double  
2715 + OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrack* object2, string variable, string function){  
2716 +  double electronMass = 0.000511;        
2717 +  double value = 0.0;    
2718 +  TLorentzVector fourVector1(0, 0, 0, 0);        
2719 +  TLorentzVector fourVector2(0, 0, 0, 0);        
2720 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));  
2721 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2722 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);    
2723 +  else if(variable == "invMass"){        
2724 +    fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, electronMass);    
2725 +    fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, electronMass );    
2726 +        
2727 +    value = (fourVector1 + fourVector2).M();    
2728 +  }
2729 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}  
2730 +  value = applyFunction(function, value);        
2731 +  return value;  
2732 + }
2733 +
2734 +
2735 +
2736 + double
2737 + OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrack* object2, string variable, string function){
2738 +  double pionMass = 0.140;
2739 +  double muonMass = 0.106;
2740 +  double value = 0.0;
2741 +  TLorentzVector fourVector1(0, 0, 0, 0);
2742 +  TLorentzVector fourVector2(0, 0, 0, 0);
2743 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2744 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2745 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2746 +  else if(variable == "invMass"){
2747 +    fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, muonMass);
2748 +    fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, pionMass );
2749 +
2750 +    value = (fourVector1 + fourVector2).M();
2751 +  }
2752 +
2753 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2754 +  value = applyFunction(function, value);
2755 +  return value;
2756 + }
2757  
2758  
2759 + double
2760 + OSUAnalysis::valueLookup (const BNtau* object1, const BNtau* object2, string variable, string function){
2761 +  double value = 0.0;
2762 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2763 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2764 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2765 +  else if(variable == "invMass"){
2766 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2767 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2768 +    value = (fourVector1 + fourVector2).M();
2769 +  }
2770 +
2771    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2772 +  value = applyFunction(function, value);
2773 +  return value;
2774 + }
2775 +
2776 +
2777 + double
2778 + OSUAnalysis::valueLookup (const BNmuon* object1, const BNtau* object2, string variable, string function){
2779 +  double value = 0.0;
2780 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2781 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
2782 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2783 +  else if(variable == "invMass"){
2784 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2785 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2786 +    value = (fourVector1 + fourVector2).M();
2787 +  }
2788  
2789 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2790    value = applyFunction(function, value);
2791 +  return value;
2792 + }
2793  
2794 + double
2795 + OSUAnalysis::valueLookup (const BNtau* object1, const BNtrack* object2, string variable, string function){
2796 +  double value = 0.0;
2797 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2798 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2799 +
2800 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2801 +  value = applyFunction(function, value);
2802    return value;
2803   }
2804  
2805  
2806 < // Calculate the number of tracks in cone of DeltaR<0.5 around track1.  
2807 < // Return true iff no other tracks are found in this cone.  
2808 < int
2806 >
2807 > double
2808 > OSUAnalysis::valueLookup (const BNtrack* object1, const BNevent* object2, string variable, string function){
2809 >
2810 >  double value = 0.0;
2811 >  double pMag = sqrt(object1->pt * object1->pt +
2812 >                     object1->pz * object1->pz);  
2813 >
2814 >  if      (variable == "numPV")                      value = object2->numPV;
2815 >  else if (variable == "caloTotDeltaRp5")            value =  (object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5);
2816 >  else if (variable == "caloTotDeltaRp5ByP")         value = ((object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5)/pMag);
2817 >  else if (variable == "caloTotDeltaRp5_RhoCorr")    value = getTrkCaloTotRhoCorr(object1);  
2818 >  else if (variable == "caloTotDeltaRp5ByP_RhoCorr") value = getTrkCaloTotRhoCorr(object1) / pMag;  
2819 >
2820 >  else { std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
2821 >
2822 >  value = applyFunction(function, value);
2823 >
2824 >  return value;
2825 >
2826 > }  
2827 >
2828 > // Calculate the number of tracks in cone of DeltaR<0.5 around track1.
2829 > // Return true iff no other tracks are found in this cone.
2830 > int
2831   OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
2832    for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
2833 <    if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.  
2834 <    double deltaRtrk = deltaR(track1->eta, track2->eta, track1->phi, track2->phi);
2833 >    if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.
2834 >    double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
2835      if(deltaRtrk < 0.5) return 0;
2836    }
2837    return 1;
2838 <  
2838 >
2839   }
2840 <  
2840 >
2841 >
2842 > double
2843 > OSUAnalysis::getTrkPtRes (const BNtrack* track1){
2844 >
2845 >  double ptTrue = getTrkPtTrue(track1, mcparticles.product());
2846 >  double PtRes = (track1->pt - ptTrue) / ptTrue;
2847 >
2848 >  return PtRes;
2849 >
2850 > }
2851 >
2852 >
2853 > double
2854 > OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
2855 >  double value = -99;
2856 >  double genDeltaRLowest = 999;
2857 >
2858 >  for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
2859 >    double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
2860 >    if (genDeltaRtemp < genDeltaRLowest) {
2861 >      genDeltaRLowest = genDeltaRtemp;
2862 >      if (genDeltaRLowest < 0.05) {   // Only consider it truth-matched if DeltaR<0.15.
2863 >        double ptTrue = genPart->pt;
2864 >        value = ptTrue;
2865 >      }
2866 >    }
2867 >  }
2868 >
2869 >  return value;
2870 >
2871 > }
2872 >
2873 > double
2874 > OSUAnalysis::getTrkCaloTotRhoCorr(const BNtrack* track) {
2875 >  // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.  
2876 >  if (!useTrackCaloRhoCorr_) return -99;  
2877 >  // if (!rhokt6CaloJetsHandle_) {
2878 >  //   cout << "ERROR [getTrkCaloTotRhoCorr]:  The collection rhokt6CaloJetsHandle is not available!" << endl;  
2879 >  //   return -99;  
2880 >  // }
2881 >  double radDeltaRCone = 0.5;  
2882 >  double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);  // Define effective area as pi*r^2, where r is radius of DeltaR cone.  
2883 >  double rawCaloTot = track->caloHadDeltaRp5 + track->caloEMDeltaRp5;  
2884 >  double caloTotRhoCorrCalo = TMath::Max(0., rawCaloTot - rhoCorr_kt6CaloJets);  
2885 >  return caloTotRhoCorrCalo;  
2886 >
2887 > }
2888 >
2889 >
2890 >
2891 >
2892   //creates a map of the dead Ecal channels in the barrel and endcap
2893   //to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
2894   //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
2895 < void
2895 > void
2896   OSUAnalysis::WriteDeadEcal (){
2897    double etaEcal, phiEcal;
2898    ifstream DeadEcalFile(deadEcalFile_);
2899    if(!DeadEcalFile) {
2900      cout << "Error: DeadEcalFile has not been found." << endl;
2901 <    return;
2902 <  }
2903 <  if(DeadEcalVec.size()!= 0){
2901 >    return;
2902 >  }
2903 >  if(DeadEcalVec.size()!= 0){
2904      cout << "Error: DeadEcalVec has a nonzero size" << endl;
2905      return;
2906 <  }
2906 >  }
2907    while(!DeadEcalFile.eof())
2908      {
2909        DeadEcalFile >> etaEcal >> phiEcal;
# Line 2356 | Line 2916 | OSUAnalysis::WriteDeadEcal (){
2916   }
2917  
2918   //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
2919 < int
2919 > int
2920   OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
2921 +  double deltaRLowest = 999;
2922    int value = 0;
2923    if (DeadEcalVec.size() == 0) WriteDeadEcal();
2924    for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
2925      double eta = ecal->etaEcal;
2926      double phi = ecal->phiEcal;
2927 <    if (deltaR(eta, track1->eta, phi, track1->phi)<0.05) {value = 1;}
2928 <    else {value = 0;}
2927 >    double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
2928 >    if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
2929    }
2930 +  if (deltaRLowest<0.05) {value = 1;}
2931 +  else {value = 0;}
2932    return value;
2933   }
2934  
2935 <
2936 <
2935 > // Returns the smallest DeltaR between the object and any generated true particle in the event.  
2936 > template <class InputObject>
2937 > double OSUAnalysis::getGenDeltaRLowest(InputObject object){
2938 >  double genDeltaRLowest = 999.;
2939 >  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2940 >    double deltaRtemp = deltaR(mcparticle->eta, mcparticle->phi, object->eta, object->phi);
2941 >    if (deltaRtemp < genDeltaRLowest) genDeltaRLowest = deltaRtemp;
2942 >  }
2943 >  return genDeltaRLowest;
2944 > }
2945  
2946   double
2947   OSUAnalysis::applyFunction(string function, double value){
2948  
2949    if(function == "abs") value = fabs(value);
2950    else if(function == "fabs") value = fabs(value);
2951 <
2951 >  else if(function == "log10") value = log10(value);
2952 >  else if(function == "log") value = log10(value);
2953  
2954    else if(function == "") value = value;
2955    else{std::cout << "WARNING: invalid function '" << function << "'\n";}
# Line 2389 | Line 2961 | OSUAnalysis::applyFunction(string functi
2961  
2962   template <class InputCollection>
2963   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2964 <
2964 >
2965 >  if (currentCut.inputCollection.find("pair")!=std::string::npos)  {
2966 >    string obj1, obj2;
2967 >    getTwoObjs(currentCut.inputCollection, obj1, obj2);
2968 >    if (inputType==obj1 ||
2969 >        inputType==obj2) {
2970 >      // Do not add a cut to individualFlags or cumulativeFlags, if the cut is on a paired collection,
2971 >      // and the inputType is a member of the pair.  
2972 >      // The cut will instead be applied when the setObjectFlags() is called for the paired collection.  
2973 >      // For example, if currentCut.inputCollection==electron-muon pairs,
2974 >      // then the flags should not be set here when inputType==muons or inputType==electrons.  
2975 >      return;
2976 >    }  
2977 >  }    
2978  
2979    for (uint object = 0; object != inputCollection->size(); object++){
2980  
2396
2981      bool decision = true;//object passes if this cut doesn't cut on that type of object
2982  
2399
2983      if(currentCut.inputCollection == inputType){
2984  
2985        vector<bool> subcutDecisions;
2986        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2987 <        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2988 <        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2987 >        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2988 >        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2989 >
2990        }
2991        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2992        else{
2993 <        bool tempDecision = true;
2994 <        for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2995 <          if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2996 <            tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2997 <          else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2993 >        bool tempDecision = true;
2994 >        for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2995 >          if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2996 >            tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2997 >          else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2998              tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
2999 <        }
3000 <        decision = tempDecision;
2999 >        }
3000 >        decision = tempDecision;
3001        }
3002      }
2419    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
3003  
3004 +    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
3005  
3006      //set flags for objects that pass each cut AND all the previous cuts
3007      bool previousCumulativeFlag = true;
# Line 2427 | Line 3011 | void OSUAnalysis::setObjectFlags(cut &cu
3011      }
3012      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
3013  
2430
3014    }
3015  
3016   }
# Line 2435 | Line 3018 | void OSUAnalysis::setObjectFlags(cut &cu
3018  
3019   template <class InputCollection1, class InputCollection2>
3020   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
3021 <                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
3021 >                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
3022  
3023  
3024    bool sameObjects = false;
3025    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
3026  
3027 +  // Get the strings for the two objects that make up the pair.  
3028 +  string obj1Type, obj2Type;
3029 +  getTwoObjs(inputType, obj1Type, obj2Type);
3030 +  bool isTwoTypesOfObject = true;
3031 +  if (obj1Type==obj2Type) isTwoTypesOfObject = false;  
3032 +
3033 +  // Initialize the flags for individual objects to all be false, if the cut is on the pair.  
3034 +  // Set them to true later, if any paired object passes (in which case both of its constituents should pass).  
3035 +  if (currentCut.inputCollection == inputType) {    
3036 +    for (uint object1 = 0; object1 != inputCollection1->size(); object1++) {
3037 +      individualFlags.at(obj1Type).at(currentCutIndex).push_back(false);
3038 +      cumulativeFlags.at(obj1Type).at(currentCutIndex).push_back(false);
3039 +    }
3040 +    if (isTwoTypesOfObject) { // Only initialize the second object if it is different from the first.  
3041 +      for (uint object2 = 0; object2 != inputCollection2->size(); object2++)  {
3042 +        individualFlags.at(obj2Type).at(currentCutIndex).push_back(false);
3043 +        cumulativeFlags.at(obj2Type).at(currentCutIndex).push_back(false);
3044 +      }
3045 +    }
3046 +  }
3047 +  
3048 +  int counter = 0;
3049  
2445  int counter = 0;  
3050    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
3051      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
3052 <      
3052 >
3053        if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
3054  
3055  
# Line 2453 | Line 3057 | void OSUAnalysis::setObjectFlags(cut &cu
3057  
3058        if(currentCut.inputCollection == inputType){
3059  
3060 <        vector<bool> subcutDecisions;
3061 <        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
3062 <          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
3063 <          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
3064 <        }
3060 >        vector<bool> subcutDecisions;
3061 >        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
3062 >          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
3063 >          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
3064 >        }
3065  
3066 <        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
3067 <        else{
3068 <          bool tempDecision = true;
3069 <          for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
3070 <            if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
3071 <              tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
3072 <            else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
3073 <              tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
3074 <          }
3075 <          decision = tempDecision;
3076 <        }
3066 >        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
3067 >        else{
3068 >          bool tempDecision = subcutDecisions.at(0);
3069 >          for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
3070 >            if(currentCut.logicalOperators.at(subcutIndex-1) == "&" || currentCut.logicalOperators.at(subcutIndex-1) == "&&")
3071 >              tempDecision = tempDecision && subcutDecisions.at(subcutIndex);
3072 >            else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
3073 >              tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
3074 >          }
3075 >          decision = tempDecision;
3076 >        }
3077        }
3078 <      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
3079 <      
3078 >      // if (decision) isPassObj1.at(object1) = true;
3079 >      // if (decision) isPassObj2.at(object2) = true;
3080 >      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);  
3081 >      if (decision && currentCut.inputCollection == inputType) {  // only set the flags for the individual objects if the pair object is being cut on  
3082 >        individualFlags.at(obj1Type).at(currentCutIndex).at(object1) = true;  
3083 >        individualFlags.at(obj2Type).at(currentCutIndex).at(object2) = true;  
3084 >      }  
3085 >
3086        //set flags for objects that pass each cut AND all the previous cuts
3087        bool previousCumulativeFlag = true;
3088        for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
3089 <        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
3090 <        else{ previousCumulativeFlag = false; break;}
3089 >        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
3090 >        else{ previousCumulativeFlag = false; break;}
3091        }
3092        //apply flags for the components of the composite object as well
3093        bool currentCumulativeFlag = true;
# Line 2486 | Line 3096 | void OSUAnalysis::setObjectFlags(cut &cu
3096        else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
3097        else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
3098        cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
3099 +      if (currentCumulativeFlag && currentCut.inputCollection == inputType) {  // only set the flags for the individual objects if the pair object is being cut on  
3100 +        cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1) = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1);  
3101 +        cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2) = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2);  
3102 +      }
3103 +      counter++;
3104  
3105 <      counter++;      
3106 <    }
2492 <    
2493 <  }
2494 <
3105 >    } // end   for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
3106 >  }  // end   for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
3107  
3108   }
3109  
3110  
3111 + bool OSUAnalysis::getPreviousCumulativeFlags(uint currentCutIndex, flagMap &individualFlags, string obj1Type, uint object1) {
3112 +  // Return true iff for the collection obj1Type, the element with index object1 has individal flags set to true for
3113 +  // all cuts up to currentCutIndex  
3114 +  bool previousCumulativeFlag = true;  
3115 +  for (uint previousCutIndex = 0; previousCutIndex < currentCutIndex; previousCutIndex++) {  
3116 +    if (previousCumulativeFlag && individualFlags.at(obj1Type).at(previousCutIndex).at(object1)) previousCumulativeFlag = true;
3117 +    else {  
3118 +      previousCumulativeFlag = false; break;  
3119 +    }
3120 +  }
3121 +  return previousCumulativeFlag;  
3122 + }  
3123 +
3124 +
3125   template <class InputCollection>
3126 < void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double puScaleFactor){
3126 > void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
3127  
3128    for (uint object = 0; object != inputCollection->size(); object++){
3129  
# Line 2507 | Line 3133 | void OSUAnalysis::fill1DHistogram(TH1* h
3133      string inputVariable = "";
3134      string function = "";
3135      if(currentString.find("(")==std::string::npos){
3136 <      inputVariable = currentString;// variable to cut on                                                                                                                                        
3136 >      inputVariable = currentString;// variable to cut on
3137      }
3138      else{
3139 <      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3140 <      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3141 <      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3139 >      function = currentString.substr(0,currentString.find("("));//function comes before the "("
3140 >      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3141 >      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3142      }
3143  
3144      double value = valueLookup(&inputCollection->at(object), inputVariable, function);
3145 <    histo->Fill(value,puScaleFactor);
3145 >    histo->Fill(value,scaleFactor);
3146  
3147 +    if (printEventInfo_) {
3148 +      // Write information about event to screen, for testing purposes.  
3149 +      cout << "  Info for event:  value for histogram " << histo->GetName() << ":  " << value << endl;  
3150 +    }
3151 +    
3152    }
2522
3153   }
3154  
3155   template <class InputCollection1, class InputCollection2>
3156 < void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double puScaleFactor){
3156 > void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
3157  
3158    bool sameObjects = false;
3159    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
3160  
3161 <  int pairCounter = 0;  
3161 >  int pairCounter = -1;
3162    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
3163      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
3164 <      
3164 >
3165        if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
3166  
3167 +      pairCounter++;
3168        //only take objects which have passed all cuts and pairs which have passed all cuts
3169        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
3170        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
# Line 2543 | Line 3174 | void OSUAnalysis::fill1DHistogram(TH1* h
3174        string inputVariable = "";
3175        string function = "";
3176        if(currentString.find("(")==std::string::npos){
3177 <        inputVariable = currentString;// variable to cut on                                                                                                          
3177 >        inputVariable = currentString;// variable to cut on
3178        }
3179        else{
3180 <        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3181 <        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3182 <        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3180 >        function = currentString.substr(0,currentString.find("("));//function comes before the "("
3181 >        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3182 >        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3183        }
3184 <      
3184 >
3185        double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
3186 <      histo->Fill(value,puScaleFactor);
3186 >      histo->Fill(value,scaleFactor);
3187  
2557      pairCounter++;      
3188      }
3189    }
3190  
# Line 2562 | Line 3192 | void OSUAnalysis::fill1DHistogram(TH1* h
3192  
3193  
3194   template <class InputCollection>
3195 < void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double puScaleFactor){
3195 > void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
3196  
3197    for (uint object = 0; object != inputCollection->size(); object++){
3198  
# Line 2572 | Line 3202 | void OSUAnalysis::fill2DHistogram(TH2* h
3202      string inputVariable = "";
3203      string function = "";
3204      if(currentString.find("(")==std::string::npos){
3205 <      inputVariable = currentString;// variable to cut on                                                                                                                                        
3205 >      inputVariable = currentString;// variable to cut on
3206      }
3207      else{
3208 <      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3209 <      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3210 <      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3208 >      function = currentString.substr(0,currentString.find("("));//function comes before the "("
3209 >      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3210 >      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3211      }
3212      double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
3213  
# Line 2585 | Line 3215 | void OSUAnalysis::fill2DHistogram(TH2* h
3215      inputVariable = "";
3216      function = "";
3217      if(currentString.find("(")==std::string::npos){
3218 <      inputVariable = currentString;// variable to cut on                                                                                                                                        
3218 >      inputVariable = currentString;// variable to cut on
3219      }
3220      else{
3221 <      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3222 <      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3223 <      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3221 >      function = currentString.substr(0,currentString.find("("));//function comes before the "("
3222 >      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3223 >      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3224      }
3225  
3226      double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
3227  
3228 <    histo->Fill(valueX,valueY,puScaleFactor);
3228 >    histo->Fill(valueX,valueY,scaleFactor);
3229  
3230    }
3231  
3232   }
3233  
3234   template <class InputCollection1, class InputCollection2>
3235 < void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double puScaleFactor){
3235 > void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
3236  
3237    bool sameObjects = false;
3238    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
3239  
3240 <  int pairCounter = 0;  
3240 >  int pairCounter = -1;
3241    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
3242      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
3243 <      
3243 >
3244        if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
3245  
3246 +      pairCounter++;
3247 +
3248        //only take objects which have passed all cuts and pairs which have passed all cuts
3249        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
3250        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
# Line 2622 | Line 3254 | void OSUAnalysis::fill2DHistogram(TH2* h
3254        string inputVariable = "";
3255        string function = "";
3256        if(currentString.find("(")==std::string::npos){
3257 <        inputVariable = currentString;// variable to cut on                                                                                                          
3257 >        inputVariable = currentString;// variable to cut on
3258        }
3259        else{
3260 <        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3261 <        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3262 <        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3260 >        function = currentString.substr(0,currentString.find("("));//function comes before the "("
3261 >        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3262 >        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3263        }
3264        double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
3265  
# Line 2635 | Line 3267 | void OSUAnalysis::fill2DHistogram(TH2* h
3267        inputVariable = "";
3268        function = "";
3269        if(currentString.find("(")==std::string::npos){
3270 <        inputVariable = currentString;// variable to cut on                                                                                                          
3270 >        inputVariable = currentString;// variable to cut on
3271        }
3272        else{
3273 <        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
3274 <        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
3275 <        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
3273 >        function = currentString.substr(0,currentString.find("("));//function comes before the "("
3274 >        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
3275 >        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
3276        }
3277        double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
3278  
3279  
3280 <      histo->Fill(valueX,valueY,puScaleFactor);
3280 >      histo->Fill(valueX,valueY,scaleFactor);
3281  
2650      pairCounter++;      
3282      }
3283    }
3284  
# Line 2664 | Line 3295 | int OSUAnalysis::getGenMatchedParticleIn
3295  
3296      double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
3297      if(currentDeltaR > 0.05) continue;
3298 <
3298 >    //     cout << std::setprecision(3) << std::setw(20)
3299 >    //          << "\tcurrentParticle:  eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
3300 >    //          << std::setw(20)
3301 >    //          << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
3302 >    //          << std::setw(20)
3303 >    //          << "\tdeltaR = " << currentDeltaR
3304 >    //          << std::setprecision(1)
3305 >    //          << std::setw(20)
3306 >    //          << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
3307 >    //          << std::setw(20)
3308 >    //          << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
3309 >    //          << std::setw(20)
3310 >    //          << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
3311      if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
3312        bestMatchIndex = mcparticle - mcparticles->begin();
3313        bestMatchDeltaR = currentDeltaR;
3314      }
3315  
3316    }
3317 <
3317 >  //   if(bestMatchDeltaR != 999)  cout << "bestMatch:  deltaR = " << bestMatchDeltaR << "   id = " << mcparticles->at(bestMatchIndex).id << "   motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
3318 >  //   else cout << "no match found..." << endl;
3319    return bestMatchIndex;
3320  
3321   }
# Line 2730 | Line 3374 | int OSUAnalysis::findTauMotherIndex(cons
3374   int OSUAnalysis::getPdgIdBinValue(int pdgId){
3375  
3376    int binValue = -999;
3377 <  int absPdgId = fabs(pdgId);  
3377 >  int absPdgId = fabs(pdgId);
3378    if(pdgId == -1) binValue = 0;
3379    else if(absPdgId <= 6 ) binValue = absPdgId;
3380    else if(absPdgId == 11 ) binValue = 7;
# Line 2756 | Line 3400 | int OSUAnalysis::getPdgIdBinValue(int pd
3400  
3401   }
3402  
3403 + const BNprimaryvertex *
3404 + OSUAnalysis::chosenVertex ()
3405 + {
3406 +  const BNprimaryvertex *chosenVertex = 0;
3407 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
3408 +    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back().size() ? cumulativeFlags.at("primaryvertexs").back() :
3409 +      cumulativeFlags.at("primaryvertexs").at(cumulativeFlags.at("primaryvertexs").size() - 2);
3410 +    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
3411 +      if(!vertexFlags.at(vertexIndex)) continue;
3412 +      chosenVertex = & primaryvertexs->at(vertexIndex);
3413 +      break;
3414 +    }
3415 +  }
3416 +  else {
3417 +    chosenVertex = &primaryvertexs->at(0);
3418 +  }
3419  
3420 < DEFINE_FWK_MODULE(OSUAnalysis);
3420 >  return chosenVertex;
3421 > }
3422 >
3423 > const BNmet *
3424 > OSUAnalysis::chosenMET ()
3425 > {
3426 >  const BNmet *chosenMET = 0;
3427 >  if(std::find(objectsToCut.begin(), objectsToCut.end(), "mets") != objectsToCut.end()) {
3428 >    vector<bool> metFlags = cumulativeFlags.at("mets").back().size() ? cumulativeFlags.at("mets").back() :
3429 >      cumulativeFlags.at("mets").at(cumulativeFlags.at("mets").size() - 2);
3430 >    for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
3431 >      if(!metFlags.at(metIndex)) continue;
3432 >      chosenMET = & mets->at(metIndex);
3433 >      break;
3434 >    }
3435 >  }
3436 >  else {
3437 >    chosenMET = &mets->at(0);
3438 >  }
3439 >
3440 >  return chosenMET;
3441 > }
3442 >
3443 > const BNelectron *
3444 > OSUAnalysis::chosenElectron ()
3445 > {
3446 >  const BNelectron *chosenElectron = 0;
3447 >  if(std::find(objectsToCut.begin(), objectsToCut.end(), "electrons") != objectsToCut.end()) {
3448 >    vector<bool> electronFlags = cumulativeFlags.at("electrons").back().size() ? cumulativeFlags.at("electrons").back() :
3449 >      cumulativeFlags.at("electrons").at(cumulativeFlags.at("electrons").size() - 2);
3450 >    for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
3451 >      if(!electronFlags.at(electronIndex)) continue;
3452 >      chosenElectron = & electrons->at(electronIndex);
3453 >      break;
3454 >    }
3455 >  }
3456 >  else {
3457 >    chosenElectron = &electrons->at(0);
3458 >  }
3459  
3460 +  return chosenElectron;
3461 + }
3462  
3463 + const BNmuon *
3464 + OSUAnalysis::chosenMuon ()
3465 + {
3466 +  const BNmuon *chosenMuon = 0;
3467 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "muons") != objectsToCut.end()) {
3468 +    vector<bool> muonFlags = cumulativeFlags.at("muons").back().size() ? cumulativeFlags.at("muons").back() :
3469 +      cumulativeFlags.at("muons").at(cumulativeFlags.at("muons").size() - 2);
3470 +    for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
3471 +      if(!muonFlags.at(muonIndex)) continue;
3472 +      chosenMuon = & muons->at(muonIndex);
3473 +      break;
3474 +    }
3475 +  }
3476 +  else {
3477 +    chosenMuon = &muons->at(0);
3478 +  }
3479 +
3480 +  return chosenMuon;
3481 + }
3482 +
3483 +
3484 + DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines