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.49 by wulsin, Wed Apr 17 15:39:43 2013 UTC vs.
Revision 1.109 by wulsin, Thu Jul 25 09:48:55 2013 UTC

# Line 1 | Line 1
1   #include "OSUT3Analysis/AnaTools/plugins/OSUAnalysis.h"
2
2   OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
3    // Retrieve parameters from the configuration file.
4    jets_ (cfg.getParameter<edm::InputTag> ("jets")),
# Line 17 | Line 16 | OSUAnalysis::OSUAnalysis (const edm::Par
16    photons_ (cfg.getParameter<edm::InputTag> ("photons")),
17    superclusters_ (cfg.getParameter<edm::InputTag> ("superclusters")),
18    triggers_ (cfg.getParameter<edm::InputTag> ("triggers")),
19 <  puFile_ (cfg.getParameter<std::string> ("puFile")),
20 <  deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")),
21 <  muonSFFile_ (cfg.getParameter<std::string> ("muonSFFile")),
22 <  dataPU_ (cfg.getParameter<std::string> ("dataPU")),
23 <  electronSFID_ (cfg.getParameter<std::string> ("electronSFID")),
24 <  muonSF_ (cfg.getParameter<std::string> ("muonSF")),
25 <  dataset_ (cfg.getParameter<std::string> ("dataset")),
26 <  datasetType_ (cfg.getParameter<std::string> ("datasetType")),
19 >  trigobjs_ (cfg.getParameter<edm::InputTag> ("trigobjs")),
20 >  puFile_ (cfg.getParameter<string> ("puFile")),
21 >  deadEcalFile_ (cfg.getParameter<string> ("deadEcalFile")),
22 >  muonSFFile_ (cfg.getParameter<string> ("muonSFFile")),
23 >  dataPU_ (cfg.getParameter<string> ("dataPU")),
24 >  electronSFID_ (cfg.getParameter<string> ("electronSFID")),
25 >  muonSF_ (cfg.getParameter<string> ("muonSF")),
26 >  dataset_ (cfg.getParameter<string> ("dataset")),
27 >  datasetType_ (cfg.getParameter<string> ("datasetType")),
28    channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
29    histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
30 +  useEDMFormat_   (cfg.getParameter<bool>("useEDMFormat")),
31 +  treeBranchSets_   (cfg.getParameter<vector<edm::ParameterSet> >("treeBranchSets")),
32    plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
33    doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting")),
34 <  printEventInfo_ (cfg.getParameter<bool> ("printEventInfo")),
34 >  applyLeptonSF_ (cfg.getParameter<bool> ("applyLeptonSF")),
35 >  applyBtagSF_ (cfg.getParameter<bool> ("applyBtagSF")),
36 >  minBtag_ (cfg.getParameter<int> ("minBtag")),
37 >  printEventInfo_      (cfg.getParameter<bool> ("printEventInfo")),
38 >  printAllTriggers_    (cfg.getParameter<bool> ("printAllTriggers")),
39    useTrackCaloRhoCorr_ (cfg.getParameter<bool> ("useTrackCaloRhoCorr")),
40 <  stopCTau_ (cfg.getParameter<vector<double> > ("stopCTau"))
40 >  stopCTau_ (cfg.getParameter<vector<double> > ("stopCTau")),
41 >  GetPlotsAfterEachCut_ (cfg.getParameter<bool> ("GetPlotsAfterEachCut")),
42 >  verbose_ (cfg.getParameter<int> ("verbose"))
43 > {
44  
45 < {
45 >  if (verbose_) printEventInfo_ = true;  
46 >  if (verbose_) clog << "Beginning OSUAnalysis::OSUAnalysis constructor." << endl;  
47  
48 <  TH1::SetDefaultSumw2 ();
48 >  TH1::SetDefaultSumw2();
49  
50    //create pile-up reweighting object, if necessary
51    if(datasetType_ != "data") {
52      if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
53 <    //    muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
54 <    //    electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
53 >    if (applyLeptonSF_){
54 >      muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
55 >      electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
56 >    }
57 >    if (applyBtagSF_){
58 >      bTagSFWeight_ = new BtagSFWeight;
59 >    }
60    }
61 < #ifdef DISPLACED_SUSY
62 <  if (datasetType_ == "signalMC")
48 <    cTauWeight_ = new CTauWeight (stopCTau_.at (0), stopCTau_.at (1), stops_);
49 < #endif
61 >  if (datasetType_ == "signalMC" && regex_match (dataset_, regex ("stop.*to.*_.*mm.*")))
62 >    stopCTauWeight_ = new StopCTauWeight (stopCTau_.at(0), stopCTau_.at(1), stops_);
63  
64  
65    // Construct Cutflow Objects. These store the results of cut decisions and
66    // handle filling cut flow histograms.
67    masterCutFlow_ = new CutFlow (fs_);
55  std::vector<TFileDirectory> directories;
68  
69    //always get vertex collection so we can assign the primary vertex in the event
70    objectsToGet.push_back("primaryvertexs");
59
60  //always make the plot of number of primary vertices (to check pile-up reweighting)
71    objectsToPlot.push_back("primaryvertexs");
72 +  objectsToFlag.push_back("primaryvertexs");
73 +
74  
75    //always get the MC particles to do GEN-matching
76    objectsToGet.push_back("mcparticles");
# Line 66 | Line 78 | OSUAnalysis::OSUAnalysis (const edm::Par
78    //always get the event collection to do pile-up reweighting
79    objectsToGet.push_back("events");
80  
81 +
82 +  // Parse the tree variable definitions.
83 +  for (uint iBranchSet = 0; !useEDMFormat_ && iBranchSet<treeBranchSets_.size(); iBranchSet++) {
84 +    string tempInputCollection = treeBranchSets_.at(iBranchSet).getParameter<string> ("inputCollection");
85 +    if(tempInputCollection.find("pairs")!=string::npos) { clog << "Warning:  tree filling is not configured for pairs of objects, so will not work for collection: " << tempInputCollection << endl; }
86 +    objectsToGet.push_back(tempInputCollection);
87 +    objectsToFlag.push_back(tempInputCollection);
88 +
89 +    vector<string> branchList(treeBranchSets_.at(iBranchSet).getParameter<vector<string> >("branches"));
90 +
91 +    for (uint iBranch = 0; iBranch<branchList.size(); iBranch++) {
92 +      BranchSpecs br;
93 +      br.inputCollection = tempInputCollection;
94 +      br.inputVariable = branchList.at(iBranch);
95 +      TString newName = TString(br.inputCollection) + "_" + TString(br.inputVariable);
96 +      br.name = string(newName.Data());
97 +      treeBranches_.push_back(br);
98 +      if (verbose_>3) clog << "   Adding branch to BNTree: " << br.name << endl;  
99 +    }
100 +
101 +  } // end   for (uint iBranchSet = 0; iBranchSet<treeBranchSets_.size(); iBranchSet++)
102 +
103 +
104    //parse the histogram definitions
105    for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
106  
107      string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
108      if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
109 +    if(tempInputCollection == "photon-muon pairs") tempInputCollection = "muon-photon pairs";
110 +    if(tempInputCollection == "photon-electron pairs") tempInputCollection = "electron-photon pairs";
111 +    if(tempInputCollection == "jet-electron pairs") tempInputCollection = "electron-jet pairs";
112 +    if(tempInputCollection == "jet-photon pairs") tempInputCollection = "photon-jet pairs";
113 +    if(tempInputCollection == "jet-muon pairs") tempInputCollection = "muon-jet pairs";
114 +    if(tempInputCollection == "event-muon pairs") tempInputCollection = "muon-event pairs";
115 +    if(tempInputCollection == "jet-met pairs")  tempInputCollection = "met-jet pairs";
116 +    if(tempInputCollection == "track-jet pairs")  tempInputCollection = "track-jet pairs";
117      if(tempInputCollection == "event-track pairs")   tempInputCollection = "track-event pairs";
118      if(tempInputCollection == "secondary muon-muon pairs")   tempInputCollection = "muon-secondary muon pairs";
119 <    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
120 <      if(tempInputCollection.find("secondary")!=std::string::npos){//secondary object
121 <        int spaceIndex = tempInputCollection.find(" ");
122 <        int secondWordLength = tempInputCollection.size() - spaceIndex;
123 <        objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
119 >    if(tempInputCollection == "secondary jet-muon pairs")   tempInputCollection = "muon-secondary jet pairs";
120 >    if(tempInputCollection == "secondary photon-muon pairs")   tempInputCollection = "muon-secondary photon pairs";
121 >    if(tempInputCollection == "secondary jet-electron pairs")   tempInputCollection = "electron-secondary jet pairs";
122 >    if(tempInputCollection == "secondary jet-photon pairs")   tempInputCollection = "photon-secondary jet pairs";
123 >    if(tempInputCollection == "secondary jet-jet pairs")   tempInputCollection = "jet-secondary jet pairs";
124 >    if(tempInputCollection == "secondary electron-electron pairs")   tempInputCollection = "electron-secondary electron pairs";
125 >    if(tempInputCollection == "trigobj-electron pairs")   tempInputCollection = "electron-trigobj pairs";
126 >    if(tempInputCollection == "trigobj-muon pairs")   tempInputCollection = "muon-trigobj pairs";
127 >    if(tempInputCollection.find("pairs")==string::npos){ //just a single object
128 >      if(tempInputCollection.find("secondary")!=string::npos){//secondary object
129 >        int spaceIndex = tempInputCollection.find(" ");
130 >        int secondWordLength = tempInputCollection.size() - spaceIndex;
131 >        objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
132        }
133        else{
134 <        objectsToGet.push_back(tempInputCollection);
134 >        objectsToGet.push_back(tempInputCollection);
135        }
136        objectsToPlot.push_back(tempInputCollection);
137 <      objectsToCut.push_back(tempInputCollection);
138 <    }
139 <    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;
137 >      objectsToFlag.push_back(tempInputCollection);
138 >    } else { //pair of objects, need to add the pair and the individual objects to the lists of things to Get/Plot/Cut
139 >      string obj1;
140        string obj2;
141        getTwoObjs(tempInputCollection, obj1, obj2);
142 <      string obj2ToGet = getObjToGet(obj2);  
143 <      objectsToCut.push_back(tempInputCollection);
144 <      objectsToCut.push_back(obj1);
145 <      objectsToCut.push_back(obj2);  
142 >      string obj2ToGet = getObjToGet(obj2);
143 >      objectsToFlag.push_back(tempInputCollection);
144 >      objectsToFlag.push_back(obj1);
145 >      objectsToFlag.push_back(obj2);
146        objectsToPlot.push_back(tempInputCollection);
147        objectsToPlot.push_back(obj1);
148 <      objectsToPlot.push_back(obj2);  
148 >      objectsToPlot.push_back(obj2);
149        objectsToGet.push_back(tempInputCollection);
150        objectsToGet.push_back(obj1);
151        objectsToGet.push_back(obj2ToGet);
152 <
153 <    }
154 <
152 >    } // end else
153 >    if (find(objectsToPlot.begin(), objectsToPlot.end(), "events") != objectsToPlot.end())
154 >      objectsToGet.push_back("jets");
155 >    
156 >    
157      vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
158  
159      for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
160  
161 +      vector<double> defaultValue;
162 +      defaultValue.push_back (-1.0);
163 +
164        histogram tempHistogram;
165        tempHistogram.inputCollection = tempInputCollection;
166        tempHistogram.name = histogramList_.at(currentHistogram).getParameter<string>("name");
167        tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
168 <      tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
168 >      tempHistogram.bins = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("bins", defaultValue);
169 >      tempHistogram.variableBinsX = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("variableBinsX", defaultValue);
170 >      tempHistogram.variableBinsY = histogramList_.at(currentHistogram).getUntrackedParameter<vector<double> >("variableBinsY", defaultValue);
171        tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
172  
173        histograms.push_back(tempHistogram);
174  
175      }
176 <  }
176 >  } //   for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++)
177 >
178    //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
179    sort( objectsToPlot.begin(), objectsToPlot.end() );
180    objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
# Line 127 | Line 185 | OSUAnalysis::OSUAnalysis (const edm::Par
185    for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
186  
187      string currentObject = objectsToPlot.at(currentObjectIndex);
188 <    if(currentObject != "muons" && currentObject != "secondary muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
188 >    if(currentObject != "muons" && currentObject != "secondary muons" && currentObject != "secondary electrons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "secondary photons"&& currentObject != "superclusters") continue;
189  
190      histogram tempIdHisto;
191      histogram tempMomIdHisto;
192      histogram tempGmaIdHisto;
193      histogram tempIdVsMomIdHisto;
194 +    histogram tempIdVsGmaIdHisto;
195  
196      tempIdHisto.inputCollection = currentObject;
197      tempMomIdHisto.inputCollection = currentObject;
198      tempGmaIdHisto.inputCollection = currentObject;
199      tempIdVsMomIdHisto.inputCollection = currentObject;
200 +    tempIdVsGmaIdHisto.inputCollection = currentObject;
201  
202      if(currentObject == "secondary muons") currentObject = "secondaryMuons";
203 +    if(currentObject == "secondary photons") currentObject = "secondaryPhotons";
204 +    if(currentObject == "secondary electrons") currentObject = "secondaryElectrons";
205  
206      currentObject = currentObject.substr(0, currentObject.size()-1);
207      tempIdHisto.name = currentObject+"GenMatchId";
208      tempMomIdHisto.name = currentObject+"GenMatchMotherId";
209      tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
210      tempIdVsMomIdHisto.name = currentObject+"GenMatchIdVsMotherId";
211 +    tempIdVsGmaIdHisto.name = currentObject+"GenMatchIdVsGrandmotherId";
212  
213      currentObject.at(0) = toupper(currentObject.at(0));
214      tempIdHisto.title = currentObject+" Gen-matched Particle";
215      tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
216      tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
217      tempIdVsMomIdHisto.title = currentObject+" Gen-matched Particle's Mother vs. Particle;Particle;Mother";
218 +    tempIdVsGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother vs. Particle;Particle;Grandmother";
219 +
220  
221 <    int maxNum = 24;
221 >    int maxNum = 25;
222      vector<double> binVector;
223      binVector.push_back(maxNum);
224      binVector.push_back(0);
# Line 171 | Line 236 | OSUAnalysis::OSUAnalysis (const edm::Par
236      tempIdVsMomIdHisto.bins = binVector;
237      tempIdVsMomIdHisto.inputVariables.push_back("genMatchedId");
238      tempIdVsMomIdHisto.inputVariables.push_back("genMatchedMotherIdReverse");
239 +    tempIdVsGmaIdHisto.bins = binVector;
240 +    tempIdVsGmaIdHisto.inputVariables.push_back("genMatchedId");
241 +    tempIdVsGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherIdReverse");
242  
243      histograms.push_back(tempIdHisto);
244      histograms.push_back(tempMomIdHisto);
245      histograms.push_back(tempGmaIdHisto);
246      histograms.push_back(tempIdVsMomIdHisto);
247 +    histograms.push_back(tempIdVsGmaIdHisto);
248    }
249  
250  
# Line 187 | Line 256 | OSUAnalysis::OSUAnalysis (const edm::Par
256      string channelName  (channels_.at(currentChannel).getParameter<string>("name"));
257      tempChannel.name = channelName;
258      TString channelLabel = channelName;
259 +    if (verbose_) clog << "Processing channel:  " << channelName << endl;  
260  
261      //set triggers for this channel
262      vector<string> triggerNames;
263      triggerNames.clear();
264 +    vector<string> triggerToVetoNames;
265 +    triggerToVetoNames.clear();
266 +
267      tempChannel.triggers.clear();
268 +    tempChannel.triggersToVeto.clear();
269      if(channels_.at(currentChannel).exists("triggers")){
270        triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
271        for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
272          tempChannel.triggers.push_back(triggerNames.at(trigger));
273        objectsToGet.push_back("triggers");
274      }
275 <
275 >    if(channels_.at(currentChannel).exists("triggersToVeto")){
276 >      triggerToVetoNames = channels_.at(currentChannel).getParameter<vector<string> >("triggersToVeto");
277 >      for(uint trigger = 0; trigger!= triggerToVetoNames.size(); trigger++)
278 >        tempChannel.triggersToVeto.push_back(triggerToVetoNames.at(trigger));
279 >      objectsToGet.push_back("triggers");
280 >    }
281  
282  
283  
284      //create cutFlow for this channel
285      cutFlows_.push_back (new CutFlow (fs_, channelName));
286 <
287 <    //book a directory in the output file with the name of the channel
288 <    TFileDirectory subDir = fs_->mkdir( channelName );
210 <    directories.push_back(subDir);
211 <
212 <    std::map<std::string, TH1D*> oneDhistoMap;
213 <    oneDHists_.push_back(oneDhistoMap);
214 <    std::map<std::string, TH2D*> twoDhistoMap;
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);
222 <      int numBinsElements = currentHistogram.bins.size();
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);
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);
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));
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));
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
270 <
271 <      vector<TString> labelArray;
272 <      labelArray.push_back("unmatched");
273 <      labelArray.push_back("u");
274 <      labelArray.push_back("d");
275 <      labelArray.push_back("s");
276 <      labelArray.push_back("c");
277 <      labelArray.push_back("b");
278 <      labelArray.push_back("t");
279 <      labelArray.push_back("e");
280 <      labelArray.push_back("#mu");
281 <      labelArray.push_back("#tau");
282 <      labelArray.push_back("#nu");
283 <      labelArray.push_back("g");
284 <      labelArray.push_back("#gamma");
285 <      labelArray.push_back("Z");
286 <      labelArray.push_back("W");
287 <      labelArray.push_back("light meson");
288 <      labelArray.push_back("K meson");
289 <      labelArray.push_back("D meson");
290 <      labelArray.push_back("B meson");
291 <      labelArray.push_back("light baryon");
292 <      labelArray.push_back("strange baryon");
293 <      labelArray.push_back("charm baryon");
294 <      labelArray.push_back("bottom baryon");
295 <      labelArray.push_back("other");
296 <
297 <      for(int bin = 0; bin !=currentHistogram.bins.at(0); 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 <
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;
320 <
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);
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);
344 <    }
345 <
346 <
347 <
348 <
286 >    vector<TFileDirectory> directories; //vector of directories in the output file.
287 >    vector<TFileDirectory> treeDirectories; //vector of directories for trees in the output file.
288 >    vector<string> subSubDirNames;//subdirectories in each channel.
289      //get list of cuts for this channel
290      vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
291  
292 +
293      //loop over and parse all cuts
294      for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
295        cut tempCut;
296        //store input collection for cut
297        string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
298        if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
299 +      if(tempInputCollection == "photon-electron pairs") tempInputCollection = "electron-photon pairs";
300 +      if(tempInputCollection == "photon-muon pairs") tempInputCollection = "muon-photon pairs";
301 +      if(tempInputCollection == "jet-electron pairs") tempInputCollection = "electron-jet pairs";
302 +      if(tempInputCollection == "jet-muon pairs") tempInputCollection = "muon-jet pairs";
303 +      if(tempInputCollection == "event-muon pairs") tempInputCollection = "muon-event pairs";
304 +      if(tempInputCollection == "jet-met pairs")  tempInputCollection = "met-jet pairs";
305 +      if(tempInputCollection == "track-jet pairs")  tempInputCollection = "track-jet pairs";
306 +      if(tempInputCollection == "jet-photon pairs") tempInputCollection = "photon-jet pairs";
307        if(tempInputCollection == "event-track pairs")   tempInputCollection = "track-event pairs";
308        if(tempInputCollection == "secondary muon-muon pairs")   tempInputCollection = "muon-secondary muon pairs";
309 +      if(tempInputCollection == "secondary jet-jet pairs")   tempInputCollection = "jet-secondary jet pairs";
310 +      if(tempInputCollection == "secondary jet-muon pairs")   tempInputCollection = "muon-secondary jet pairs";
311 +      if(tempInputCollection == "secondary photon-muon pairs")   tempInputCollection = "muon-secondary photon pairs";
312 +      if(tempInputCollection == "secondary jet-photon pairs")   tempInputCollection = "photon-secondary jet pairs";
313 +      if(tempInputCollection == "secondary jet-electron pairs")   tempInputCollection = "electron-secondary jet pairs";
314 +      if(tempInputCollection == "secondary electron-electron pairs")   tempInputCollection = "electron-secondary electron pairs";
315 +      if(tempInputCollection == "trigobj-electron pairs")   tempInputCollection = "electron-trigobj pairs";
316 +      if(tempInputCollection == "trigobj-muon pairs")   tempInputCollection = "muon-trigobj pairs";
317        tempCut.inputCollection = tempInputCollection;
318 <      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
319 <        if(tempInputCollection.find("secondary")!=std::string::npos){//secondary object
320 <          int spaceIndex = tempInputCollection.find(" ");
321 <          int secondWordLength = tempInputCollection.size() - spaceIndex;
322 <          objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
323 <        }
324 <        else{
325 <          objectsToGet.push_back(tempInputCollection);
326 <        }
327 <        objectsToCut.push_back(tempInputCollection);
318 >      if(tempInputCollection.find("pairs")==string::npos){ //just a single object
319 >        if(tempInputCollection.find("secondary")!=string::npos){//secondary object
320 >          int spaceIndex = tempInputCollection.find(" ");
321 >          int secondWordLength = tempInputCollection.size() - spaceIndex;
322 >          objectsToGet.push_back(tempInputCollection.substr(spaceIndex+1,secondWordLength));
323 >        }
324 >        else{
325 >          objectsToGet.push_back(tempInputCollection);
326 >        }
327 >        objectsToCut.push_back(tempInputCollection);
328 >        objectsToFlag.push_back(tempInputCollection);
329        }
330        else{//pair of objects, need to add them both to objectsToGet
331 <        string obj1;
332 <        string obj2;
333 <        getTwoObjs(tempInputCollection, obj1, obj2);
334 <        string obj2ToGet = getObjToGet(obj2);  
331 >        string obj1;
332 >        string obj2;
333 >        getTwoObjs(tempInputCollection, obj1, obj2);
334 >        string obj2ToGet = getObjToGet(obj2);
335          objectsToCut.push_back(tempInputCollection);
336          objectsToCut.push_back(obj1);
337 <        objectsToCut.push_back(obj2);  
337 >        objectsToCut.push_back(obj2);
338 >        objectsToFlag.push_back(tempInputCollection);
339 >        objectsToFlag.push_back(obj1);
340 >        objectsToFlag.push_back(obj2);
341          objectsToGet.push_back(tempInputCollection);
342          objectsToGet.push_back(obj1);
343          objectsToGet.push_back(obj2ToGet);
# Line 387 | Line 348 | OSUAnalysis::OSUAnalysis (const edm::Par
348  
349        //split cut string into parts and store them
350        string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
351 <      std::vector<string> cutStringVector = splitString(cutString);
351 >      vector<string> cutStringVector = splitString(cutString);
352        if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
353 <        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
353 >        clog << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
354          exit(0);
355        }
356        tempCut.numSubcuts = (cutStringVector.size()+1)/4;
357        for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
358          int indexOffset = 4 * subcutIndex;
359          string currentVariableString = cutStringVector.at(indexOffset);
360 <        if(currentVariableString.find("(")==std::string::npos){
360 >        if(currentVariableString.find("(")==string::npos){
361            tempCut.functions.push_back("");//no function was specified
362            tempCut.variables.push_back(currentVariableString);// variable to cut on
363          }
# Line 407 | Line 368 | OSUAnalysis::OSUAnalysis (const edm::Par
368          }
369          tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
370          tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
371 +        tempCut.cutStringValues.push_back(cutStringVector.at(indexOffset+2));// string value to pass cut
372          if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
373        }
374  
375        //get number of objects required to pass cut for event to pass
376        string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
377 <      std::vector<string> numberRequiredVector = splitString(numberRequiredString);
377 >      vector<string> numberRequiredVector = splitString(numberRequiredString);
378        if(numberRequiredVector.size()!=2){
379 <        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
379 >        clog << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
380          exit(0);
381        }
382  
# Line 422 | Line 384 | OSUAnalysis::OSUAnalysis (const edm::Par
384        tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
385        tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
386  
387 <
387 >      //Set up vectors to store the directories and subDirectories for each channel.
388 >      string subSubDirName;
389        string tempCutName;
390        if(cuts_.at(currentCut).exists("alias")){
391          tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
392 +        subSubDirName = "After " + tempCutName + " Cut Applied";
393        }
394        else{
395          //construct string for cutflow table
# Line 433 | Line 397 | OSUAnalysis::OSUAnalysis (const edm::Par
397          string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
398          string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
399          tempCutName = cutName;
400 +        subSubDirName = "After " + numberRequiredString + " " + collectionString + " with " + cutString + " Cut Applied";
401 +      }
402 +
403 +      tempCut.isVeto = false;
404 +      if(cuts_.at(currentCut).exists("isVeto")){
405 +        bool isVeto = cuts_.at(currentCut).getParameter<bool> ("isVeto");
406 +        if (isVeto == true) tempCut.isVeto = true;
407        }
408 +      subSubDirNames.push_back(subSubDirName);
409        tempCut.name = tempCutName;
410 +      if (verbose_) clog << "Setting up cut, index: " << currentCut << ", input collection: " << tempInputCollection<< ", name: " << tempCut.name << endl;  
411  
412        tempChannel.cuts.push_back(tempCut);
413  
414  
415      }//end loop over cuts
416  
417 +    vector<map<string, TH1D*>> oneDHistsTmp;
418 +    vector<map<string, TH2D*>> twoDHistsTmp;
419 +    //book a directory in the output file with the name of the channel
420 +    TFileDirectory subDir = fs_->mkdir( channelName );
421 +    //loop over the cuts to set up subdirectory for each cut if GetPlotsAfterEachCut_ is true.
422 +    if(GetPlotsAfterEachCut_){
423 +      for( uint currentDir=0; currentDir != subSubDirNames.size(); currentDir++){
424 +        TFileDirectory subSubDir = subDir.mkdir( subSubDirNames[currentDir] );
425 +        directories.push_back(subSubDir);
426 +        map<string, TH1D*> oneDhistoMap;
427 +        oneDHistsTmp.push_back(oneDhistoMap);
428 +        map<string, TH2D*> twoDhistoMap;
429 +        twoDHistsTmp.push_back(twoDhistoMap);
430 +      }
431 +      treeDirectories.push_back(subDir);
432 +      oneDHists_.push_back(oneDHistsTmp);
433 +      twoDHists_.push_back(twoDHistsTmp);
434 +    }
435 +    //only set up directories with names of the channels if GetPlotsAfterEachCut_ is false.
436 +    else{
437 +      map<string, TH1D*> oneDhistoMap;
438 +      oneDHistsTmp.push_back(oneDhistoMap);
439 +      map<string, TH2D*> twoDhistoMap;
440 +      twoDHistsTmp.push_back(twoDhistoMap);
441 +      oneDHists_.push_back(oneDHistsTmp);
442 +      twoDHists_.push_back(twoDHistsTmp);
443 +      directories.push_back(subDir);
444 +      treeDirectories.push_back(subDir);
445 +
446 +    }
447 +
448 +    for(uint currentDir = 0; !useEDMFormat_ && currentDir != treeDirectories.size(); currentDir++){
449 +      TTree* newTree = treeDirectories.at(currentDir).make<TTree> (TString("BNTree_"+channelLabel), TString("BNTree_"+channelLabel));
450 +      BNTrees_.push_back(newTree);
451 +      for (uint iBranch = 0; iBranch < treeBranches_.size(); iBranch++){
452 +        BranchSpecs currentVar = treeBranches_.at(iBranch);
453 +        vector<float> newVec;
454 +        BNTreeBranchVals_[currentVar.name] = newVec;
455 +        BNTrees_.back()->Branch(TString(currentVar.name), &BNTreeBranchVals_.at(currentVar.name));
456 +      } // end for (uint iBranch = 0; iBranch < treeBranches_.size(); iBranch++)
457 +    }
458 +
459 +    //book all histograms included in the configuration
460 +    for(uint currentDir = 0; currentDir != directories.size(); currentDir++){//loop over all the directories.
461 +
462 +      for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
463 +
464 +        histogram currentHistogram = histograms.at(currentHistogramIndex);
465 +        int numBinsElements = currentHistogram.bins.size();
466 +        int numInputVariables = currentHistogram.inputVariables.size();
467 +        int numBinEdgesX = currentHistogram.variableBinsX.size();
468 +        int numBinEdgesY = currentHistogram.variableBinsY.size();
469 +
470 +        if(numBinsElements == 1){
471 +          if(numBinEdgesX > 1){
472 +            if(numBinEdgesY > 1)
473 +              numBinsElements = 6;
474 +            else
475 +              numBinsElements = 3;
476 +          }
477 +        }
478 +        if(numBinsElements != 3 && numBinsElements !=6){
479 +          clog << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
480 +          exit(0);
481 +        }
482 +        else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
483 +          clog << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
484 +          exit(0);
485 +        }
486 +        else if(numBinsElements == 3){
487 +          if (currentHistogram.bins.size () == 3)
488 +            oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH1D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, currentHistogram.bins.at(0), currentHistogram.bins.at(1), currentHistogram.bins.at(2));
489 +          else
490 +            {
491 +              oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH1D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, numBinEdgesX - 1, currentHistogram.variableBinsX.data ());
492 +            }
493 +        }
494 +        else if(numBinsElements == 6){
495 +          if (currentHistogram.bins.size () == 6)
496 +            twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH2D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, currentHistogram.bins.at(0), currentHistogram.bins.at(1), currentHistogram.bins.at(2),currentHistogram.bins.at(3),currentHistogram.bins.at(4),currentHistogram.bins.at(5));
497 +          else
498 +            twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name] = directories.at(currentDir).make<TH2D> (TString(currentHistogram.name),channelLabel+" channel: "+currentHistogram.title, numBinEdgesX - 1, currentHistogram.variableBinsX.data (), numBinEdgesY - 1, currentHistogram.variableBinsY.data ());
499 +        }
500 +
501 +
502 +        if(currentHistogram.name.find("GenMatch")==string::npos) continue;
503 +
504 +        // bin      particle type
505 +        // ---      -------------
506 +        //  0        unmatched
507 +        //  1        u
508 +        //  2        d
509 +        //  3        s
510 +        //  4        c
511 +        //  5        b
512 +        //  6        t
513 +        //  7        e
514 +        //  8        mu
515 +        //  9        tau
516 +        // 10        nu
517 +        // 11        g
518 +        // 12        gamma
519 +        // 13        Z
520 +        // 14        W
521 +        // 15        light meson
522 +        // 16        K meson
523 +        // 17        D meson
524 +        // 18        B meson
525 +        // 19        light baryon
526 +        // 20        strange baryon
527 +        // 21        charm baryon
528 +        // 22        bottom baryon
529 +        // 23        QCD string
530 +        // 24        other
531 +
532 +        vector<TString> labelArray;
533 +        labelArray.push_back("unmatched");
534 +        labelArray.push_back("u");
535 +        labelArray.push_back("d");
536 +        labelArray.push_back("s");
537 +        labelArray.push_back("c");
538 +        labelArray.push_back("b");
539 +        labelArray.push_back("t");
540 +        labelArray.push_back("e");
541 +        labelArray.push_back("#mu");
542 +        labelArray.push_back("#tau");
543 +        labelArray.push_back("#nu");
544 +        labelArray.push_back("g");
545 +        labelArray.push_back("#gamma");
546 +        labelArray.push_back("Z");
547 +        labelArray.push_back("W");
548 +        labelArray.push_back("light meson");
549 +        labelArray.push_back("K meson");
550 +        labelArray.push_back("D meson");
551 +        labelArray.push_back("B meson");
552 +        labelArray.push_back("light baryon");
553 +        labelArray.push_back("strange baryon");
554 +        labelArray.push_back("charm baryon");
555 +        labelArray.push_back("bottom baryon");
556 +        labelArray.push_back("QCD string");
557 +        labelArray.push_back("other");
558 +
559 +        for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
560 +          if(currentHistogram.name.find("GenMatchIdVsMotherId")==string::npos && currentHistogram.name.find("GenMatchIdVsGrandmotherId")==string::npos) {
561 +            oneDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
562 +          }
563 +          else {
564 +            twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetYaxis()->SetBinLabel(bin+1,labelArray.at(currentHistogram.bins.at(0)-bin-1));
565 +            twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
566 +          }
567 +        }
568 +        if(currentHistogram.name.find("GenMatchIdVsMotherId")!=string::npos || currentHistogram.name.find("GenMatchIdVsGrandmotherId")!=string::npos) {
569 +          twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetXaxis()->CenterTitle();
570 +          twoDHists_.at(currentChannel).at(currentDir)[currentHistogram.name]->GetYaxis()->CenterTitle();
571 +        }
572 +
573 +      }  // end      for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++)
574 +
575 +
576 +      // Book a histogram for the number of each object type to be plotted.
577 +      // Name of objectToPlot here must match the name specified in OSUAnalysis::analyze().
578 +      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
579 +        string currentObject = objectsToPlot.at(currentObjectIndex);
580 +        int maxNum = 10;
581 +        if(currentObject == "mcparticles") maxNum = 50;
582 +        else if(currentObject == "primaryvertexs") maxNum = 50;
583 +
584 +        if(currentObject == "muon-muon pairs")                currentObject = "dimuonPairs";
585 +        else if(currentObject == "electron-electron pairs")   currentObject = "dielectronPairs";
586 +        else if(currentObject == "electron-muon pairs")       currentObject = "electronMuonPairs";
587 +        else if(currentObject == "electron-photon pairs")     currentObject = "electronPhotonPairs";
588 +        else if(currentObject == "muon-photon pairs")         currentObject = "muonPhotonPairs";
589 +        else if(currentObject == "secondary jets")            currentObject = "secondaryJets";
590 +        else if(currentObject == "secondary photons")            currentObject = "secondaryPhotons";
591 +        else if(currentObject == "jet-jet pairs")             currentObject = "dijetPairs";
592 +        else if(currentObject == "jet-secondary jet pairs")   currentObject = "jetSecondaryJetPairs";
593 +        else if(currentObject == "electron-jet pairs")        currentObject = "electronJetPairs";
594 +        else if(currentObject == "muon-jet pairs")            currentObject = "muonJetPairs";
595 +        else if(currentObject == "muon-event pairs")            currentObject = "muonEventPairs";
596 +        else if(currentObject == "photon-jet pairs")            currentObject = "photonJetPairs";
597 +        else if(currentObject == "met-jet pairs")             currentObject = "metJetPairs";
598 +        else if(currentObject == "track-jet pairs")           currentObject = "trackJetPairs";
599 +        else if(currentObject == "muon-secondary jet pairs")  currentObject = "muonSecondaryJetPairs";
600 +        else if(currentObject == "muon-secondary photon pairs")  currentObject = "muonSecondaryPhotonPairs";
601 +        else if(currentObject == "photon-secondary jet pairs")  currentObject = "photonSecondaryJetPairs";
602 +        else if(currentObject == "electron-secondary jet pairs")  currentObject = "electronSecondaryJetPairs";
603 +        else if(currentObject == "track-event pairs")         currentObject = "trackEventPairs";
604 +        else if(currentObject == "electron-track pairs")      currentObject = "electronTrackPairs";
605 +        else if(currentObject == "muon-track pairs")          currentObject = "muonTrackPairs";
606 +        else if(currentObject == "muon-tau pairs")            currentObject = "muonTauPairs";
607 +        else if(currentObject == "tau-tau pairs")             currentObject = "ditauPairs";
608 +        else if(currentObject == "tau-track pairs")           currentObject = "tauTrackPairs";
609 +        else if(currentObject == "muon-secondary muon pairs") currentObject = "muonSecondaryMuonPairs";
610 +        else if(currentObject == "secondary muons")           currentObject = "secondaryMuons";
611 +        else if(currentObject == "electron-secondary electron pairs") currentObject = "electronSecondaryElectronPairs";
612 +        else if(currentObject == "secondary electrons")       currentObject = "secondaryElectrons";
613 +        else if(currentObject == "electron-trigobj pairs")    currentObject = "electronTrigobjPairs";
614 +        else if(currentObject == "muon-trigobj pairs")        currentObject = "muonTrigobjPairs";
615 +
616 +        currentObject.at(0) = toupper(currentObject.at(0));
617 +        string histoName = "num" + currentObject;
618 +
619 +        if(histoName == "numPrimaryvertexs"){
620 +          string newHistoName = histoName + "BeforePileupCorrection";
621 +          oneDHists_.at(currentChannel).at(currentDir)[newHistoName] = directories.at(currentDir).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
622 +          newHistoName = histoName + "AfterPileupCorrection";
623 +          oneDHists_.at(currentChannel).at(currentDir)[newHistoName] = directories.at(currentDir).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
624 +        }
625 +        else
626 +          oneDHists_.at(currentChannel).at(currentDir)[histoName] = directories.at(currentDir).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
627 +      }
628 +
629 +      if (verbose_) {
630 +        clog << "List of 1D histograms booked for channel " << channelName << " and directory " << currentDir << ":" << endl;
631 +        for(map<string, TH1D*>::iterator
632 +              iter = oneDHists_.at(currentChannel).at(currentDir).begin();
633 +            iter  != oneDHists_.at(currentChannel).at(currentDir).end();
634 +            iter++) {
635 +          clog << " " << iter->first
636 +               << ":  name=" << iter->second->GetName()
637 +               << ", title=" << iter->second->GetTitle()
638 +               << ", bins=(" << iter->second->GetNbinsX()
639 +               << "," << iter->second->GetXaxis()->GetXmin()
640 +               << "," << iter->second->GetXaxis()->GetXmax()
641 +               << ")" << endl;
642 +        }  
643 +        clog << "List of 2D histograms booked for channel " << channelName << " and directory " << currentDir << ":" << endl;
644 +        for(map<string, TH2D*>::iterator
645 +              iter = twoDHists_.at(currentChannel).at(currentDir).begin();
646 +            iter  != twoDHists_.at(currentChannel).at(currentDir).end();
647 +            iter++) {
648 +          clog << " " << iter->first
649 +               << ":  name=" << iter->second->GetName()
650 +               << ", title=" << iter->second->GetTitle()
651 +               << ", binsX=(" << iter->second->GetNbinsX()
652 +               << "," << iter->second->GetXaxis()->GetXmin()
653 +               << "," << iter->second->GetXaxis()->GetXmax()
654 +               << ")"
655 +               << ", binsY=(" << iter->second->GetNbinsY()
656 +               << "," << iter->second->GetYaxis()->GetXmin()
657 +               << "," << iter->second->GetYaxis()->GetXmax()
658 +               << ")"
659 +               << endl;
660 +        }  
661 +      }
662 +
663 +    }//end of loop over directories
664      channels.push_back(tempChannel);
665      tempChannel.cuts.clear();
446
666    }//end loop over channels
667  
668  
669 +
670    //make unique vector of objects we need to get from the event
671    sort( objectsToGet.begin(), objectsToGet.end() );
672    objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
673    //make unique vector of objects we need to cut on
674    sort( objectsToCut.begin(), objectsToCut.end() );
675    objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
676 +  //make unique vector of objects we need to set flags for
677 +  sort( objectsToFlag.begin(), objectsToFlag.end() );
678 +  objectsToFlag.erase( unique( objectsToFlag.begin(), objectsToFlag.end() ), objectsToFlag.end() );
679  
680 +  // Move all paired objects to the end of the list, so that the flags for single objects are always set before those of paired objects.  
681 +  for (uint i=0; i<objectsToFlag.size(); i++) {
682 +    if (objectsToFlag.at(i).find("pairs")!=string::npos) { // if it contains "pairs"
683 +      objectsToFlag.push_back(objectsToFlag.at(i));  
684 +      objectsToFlag.erase(objectsToFlag.begin()+i);  
685 +    }
686 +  }
687 +
688 +  if (verbose_) {
689 +    clog << "List of channels:" << endl;
690 +    for (uint i=0; i<channels.size(); i++) {
691 +      clog << " " << channels.at(i).name << endl;
692 +    }
693 +    clog << "List of objects to get:" << endl;
694 +    for (uint i=0; i<objectsToGet.size(); i++) {
695 +      clog << " " << objectsToGet.at(i) << endl;
696 +    }
697 +    clog << "List of objects to plot:" << endl;
698 +    for (uint i=0; i<objectsToPlot.size(); i++) {
699 +      clog << " " << objectsToPlot.at(i) << endl;
700 +    }
701 +    clog << "List of objects to cut:" << endl;
702 +    for (uint i=0; i<objectsToCut.size(); i++) {
703 +      clog << " " << objectsToCut.at(i) << endl;
704 +    }
705 +    clog << "List of objects to flag:" << endl;
706 +    for (uint i=0; i<objectsToFlag.size(); i++) {
707 +      clog << " " << objectsToFlag.at(i) << endl;
708 +    }
709 +  }
710 +
711 +  produces<map<string, bool> > ("channelDecisions");
712 +
713 +  if (printEventInfo_) {
714 +    findEventsLog = new ofstream();  
715 +    findEventsLog->open("findEvents.txt");  
716 +    clog << "Listing run:lumi:event in file findEvents.txt for events that pass full selection (of any channel)." << endl;  
717 +  } else {
718 +    findEventsLog = 0;
719 +  }  
720 +
721 +  isFirstEvent_ = true;  
722 +
723 +  if (verbose_) clog << "Finished OSUAnalysis::OSUAnalysis constructor." << endl;  
724 +
725 +
726 +  } // end constructor OSUAnalysis::OSUAnalysis()
727  
458 }
728  
729   OSUAnalysis::~OSUAnalysis ()
730   {
731 +
732 +  if (verbose_) clog << "Beginning OSUAnalysis::OSUAnalysis destructor." << endl;  
733 +
734    // Destroying the CutFlow objects causes the cut flow numbers and time
735    // information to be printed to standard output.
736    for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
737      delete cutFlows_.at(currentChannel);
738    }
739 +
740 +  if (printEventInfo_ && findEventsLog) findEventsLog->close();  
741 +
742 +  clog << "=============================================" << endl;  
743 +  clog << "Successfully completed OSUAnalysis." << endl;  
744 +  clog << "=============================================" << endl;  
745 +
746 +  if (verbose_) clog << "Finished OSUAnalysis::OSUAnalysis destructor." << endl;  
747 +
748   }
749  
750   void
751 < OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
751 > OSUAnalysis::produce (edm::Event &event, const edm::EventSetup &setup)
752   {
472
753    // Retrieve necessary collections from the event.
754  
755 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
755 >  if (verbose_) clog << "Beginning OSUAnalysis::produce." << endl;  
756 >
757 >  if (find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
758      event.getByLabel (triggers_, triggers);
759  
760 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
760 >  if (find(objectsToGet.begin(), objectsToGet.end(), "trigobjs") != objectsToGet.end())
761 >    event.getByLabel (trigobjs_, trigobjs);
762 >
763 >  if (find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
764      event.getByLabel (jets_, jets);
765  
766 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
766 >  if (find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
767      event.getByLabel (muons_, muons);
768  
769 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
769 >  if (find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
770      event.getByLabel (electrons_, electrons);
771  
772 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
772 >  if (find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
773      event.getByLabel (events_, events);
774  
775 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
775 >  if (find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
776      event.getByLabel (taus_, taus);
777  
778 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
778 >  if (find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
779      event.getByLabel (mets_, mets);
780  
781 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
781 >  if (find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
782      event.getByLabel (tracks_, tracks);
783  
784 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
784 >  if (find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
785      event.getByLabel (genjets_, genjets);
786  
787 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
787 >  if (find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
788      event.getByLabel (mcparticles_, mcparticles);
789  
790 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
790 >  if (find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
791      event.getByLabel (primaryvertexs_, primaryvertexs);
792  
793 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
793 >  if (find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
794      event.getByLabel (bxlumis_, bxlumis);
795  
796 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
796 >  if (find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
797      event.getByLabel (photons_, photons);
798  
799 <  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
799 >  if (find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
800      event.getByLabel (superclusters_, superclusters);
801  
802 <  if (useTrackCaloRhoCorr_) {  
803 <    // Used only for pile-up correction of by-hand calculation of isolation energy.  
804 <    // This rho collection is not available in all BEANs.  
802 >  if (datasetType_ == "signalMC"){
803 >    if (find(objectsToGet.begin(), objectsToGet.end(), "stops") != objectsToGet.end())
804 >      event.getByLabel (stops_, stops);
805 >  }
806 >
807 >  if (useTrackCaloRhoCorr_) {
808 >    // Used only for pile-up correction of by-hand calculation of isolation energy.
809 >    // This rho collection is not available in all BEANs.
810      // For description of rho values for different jet reconstruction algorithms, see
811 <    // https://twiki.cern.ch/twiki/bin/view/CMS/JetAlgorithms#Algorithms  
812 <    event.getByLabel ("kt6CaloJets","rho", rhokt6CaloJetsHandle_);
811 >    // https://twiki.cern.ch/twiki/bin/view/CMS/JetAlgorithms#Algorithms
812 >    event.getByLabel ("kt6CaloJets","rho", rhokt6CaloJetsHandle_);
813    }
814  
815 +  double masterScaleFactor = 1.0;
816  
817    //get pile-up event weight
818 <  double scaleFactor = 1.00;
819 <  if(doPileupReweighting_ && datasetType_ != "data")
820 <    scaleFactor = puWeight_->at (events->at (0).numTruePV);
821 <
822 <  cTauScaleFactor_ = 1.0;
823 < #ifdef DISPLACED_SUSY
824 <  if (datasetType_ == "signalMC")
825 <    cTauScaleFactor_ = cTauWeight_->at (event);
826 < #endif
827 <  scaleFactor *= cTauScaleFactor_;
818 >  if (doPileupReweighting_ && datasetType_ != "data") {
819 >    //for "data" datasets, the numTruePV is always set to -1
820 >    if (events->at(0).numTruePV < 0 && isFirstEvent_) clog << "WARNING[OSUAnalysis::analyze]: Event has numTruePV<0.  Turning off pile-up reweighting." << endl;
821 >    else masterScaleFactor *= puWeight_->at (events->at (0).numTruePV);
822 >  }
823 >
824 >  stopCTauScaleFactor_ = 1.0;
825 >  if (datasetType_ == "signalMC" && regex_match (dataset_, regex ("stop.*to.*_.*mm.*")))
826 >    stopCTauScaleFactor_ = stopCTauWeight_->at (event);
827 >  masterScaleFactor *= stopCTauScaleFactor_;
828  
829    //loop over all channels
830  
831 +  auto_ptr<map<string, bool> > channelDecisions (new map<string, bool>);
832    for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
833      channel currentChannel = channels.at(currentChannelIndex);
834 +    if (verbose_>1) clog << " Processing channel " << currentChannel.name << endl;
835  
836      flagMap individualFlags;
837      counterMap passingCounter;
838      cumulativeFlags.clear ();
839  
840 +    for (map<string, vector<float>>::iterator iter = BNTreeBranchVals_.begin();
841 +         iter != BNTreeBranchVals_.end(); iter++) {
842 +      iter->second.clear();  // clear array
843 +    }
844 +
845      bool triggerDecision = true;
846 <    if(currentChannel.triggers.size() != 0){  //triggers specified
847 <      triggerDecision = evaluateTriggers(currentChannel.triggers,triggers.product());
846 >    if(currentChannel.triggers.size() != 0 || currentChannel.triggersToVeto.size() != 0){  //triggers specified
847 >      triggerDecision = evaluateTriggers(currentChannel.triggers, currentChannel.triggersToVeto, triggers.product());
848        cutFlows_.at(currentChannelIndex)->at ("trigger") = triggerDecision;
849      }
850  
851      //loop over all cuts
554
555
556
852      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
853        cut currentCut = currentChannel.cuts.at(currentCutIndex);
854 +      if (verbose_>2) clog << "  Processing cut, index: " << currentCutIndex << ", input collection: " << currentCut.inputCollection << ", name: " << currentCut.name << endl;  
855  
856 <      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
857 <        string currentObject = objectsToCut.at(currentObjectIndex);
856 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToFlag.size(); currentObjectIndex++){
857 >        string currentObject = objectsToFlag.at(currentObjectIndex);
858  
859          //initialize maps to get ready to set cuts
860 <        individualFlags[currentObject].push_back (vector<bool> ());
861 <        cumulativeFlags[currentObject].push_back (vector<bool> ());
860 >        individualFlags[currentObject].push_back (vector<pair<bool,bool> > ());
861 >        cumulativeFlags[currentObject].push_back (vector<pair<bool,bool> > ());
862  
863        }
568      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
569        string currentObject = objectsToCut.at(currentObjectIndex);
864  
865 <        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
865 >      //set flags for all relevant objects
866 >      for(int currentObjectIndex = -1; currentObjectIndex != int(objectsToFlag.size()); currentObjectIndex++){
867 >        string currentObject;  
868 >        if (currentObjectIndex < 0) currentObject = currentCut.inputCollection;  // In the first loop, set the flags for the collection that the cut is acting on.  That way other paired collections can access the correct flag for the current cut.  
869 >        else currentObject = objectsToFlag.at(currentObjectIndex);
870 >        if (currentObjectIndex >= 0 && currentObject == currentCut.inputCollection) continue;  // Flags have already been set for the inputCollection object, so should not be set again.  
871 >
872 >        // single object collections
873 >        if     (currentObject == "jets")                setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
874 >        else if(currentObject == "secondary jets")      setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"secondary jets");
875 >        else if(currentObject == "secondary photons")   setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"secondary photons");
876 >        else if(currentObject == "muons")               setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
877 >        else if(currentObject == "secondary muons")     setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"secondary muons");
878 >        else if(currentObject == "secondary electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"secondary electrons");
879 >        else if(currentObject == "electrons")           setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
880 >        else if(currentObject == "events")              setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
881 >        else if(currentObject == "taus")                setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
882 >        else if(currentObject == "mets")                setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
883 >        else if(currentObject == "tracks")              setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
884 >        else if(currentObject == "genjets")             setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
885 >        else if(currentObject == "mcparticles")         setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
886 >        else if(currentObject == "primaryvertexs")      setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
887 >        else if(currentObject == "bxlumis")             setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
888 >        else if(currentObject == "photons")             setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
889 >        else if(currentObject == "superclusters")       setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
890 >        else if(currentObject == "trigobjs")            setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,trigobjs.product(),"trigobjs");
891 >
892 >
893 >        // paired object collections  
894 >        else if(currentObject == "muon-muon pairs")                   setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), "muon-muon pairs");
895 >        else if(currentObject == "muon-secondary muon pairs")         setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), "muon-secondary muon pairs");
896 >
897 >        else if(currentObject == "muon-secondary photon pairs")       setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),photons.product(), "muon-secondary photon pairs");
898 >        else if(currentObject == "muon-secondary jet pairs")          setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),jets.product(), "muon-secondary jet pairs");
899 >        else if(currentObject == "photon-secondary jet pairs")        setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),jets.product(), "photon-secondary jet pairs");
900 >        else if(currentObject == "electron-secondary jet pairs")      setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),jets.product(), "electron-secondary jet pairs");
901 >        else if(currentObject == "electron-secondary electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), "electron-secondary electron pairs");
902 >
903 >        else if(currentObject == "electron-electron pairs")           setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), "electron-electron pairs");
904 >        else if(currentObject == "electron-muon pairs")     setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), "electron-muon pairs");
905 >        else if(currentObject == "jet-jet pairs")           setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),jets.product(), "jet-jet pairs");
906 >        else if(currentObject == "jet-secondary jet pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),jets.product(), "jet-secondary jet pairs");
907 >        else if(currentObject == "electron-jet pairs")      setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),jets.product(), "electron-jet pairs");
908 >        else if(currentObject == "electron-photon pairs")   setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),photons.product(), "electron-photon pairs");
909 >        else if(currentObject == "photon-jet pairs")        setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),jets.product(), "photon-jet pairs");
910 >        else if(currentObject == "muon-jet pairs")          setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),jets.product(), "muon-jet pairs");
911 >        else if(currentObject == "muon-event pairs")        setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),events.product(), "muon-event pairs");
912 >        else if(currentObject == "met-jet pairs")           setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),jets.product(), "met-jet pairs");
913 >        else if(currentObject == "track-jet pairs")         setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),jets.product(), "track-jet pairs");
914 >        else if(currentObject == "muon-photon pairs")       setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),photons.product(), "muon-photon pairs");
915 >        else if(currentObject == "track-event pairs")       setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),events.product(), "track-event pairs");
916 >        else if(currentObject == "electron-track pairs")    setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),tracks.product(),"electron-track pairs");
917 >        else if(currentObject == "muon-track pairs")        setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),tracks.product(),"muon-track pairs");
918 >        else if(currentObject == "muon-tau pairs")          setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),taus.product(),"muon-tau pairs");
919 >        else if(currentObject == "tau-tau pairs")           setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),taus.product(),"tau-tau pairs");
920 >        else if(currentObject == "tau-track pairs")         setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus .product(),tracks.product(),"tau-track pairs");
921 >        else if(currentObject == "electron-trigobj pairs")  setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),trigobjs.product(),"electron-trigobj pairs");
922 >        else if(currentObject == "muon-trigobj pairs")      setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),trigobjs.product(),"muon-trigobj pairs");
923  
924 <
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");
582 <        else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
583 <        else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
584 <        else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
585 <        else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
586 <        else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
587 <        else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
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 <
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 <        
924 >        if(currentObject == "stops" && datasetType_ == "signalMC") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,stops.product(),"stops");
925        }
926  
638
639
927      }//end loop over all cuts
928  
929  
643    //use cumulative flags to apply cuts at event level
930  
931 +    //use cumulative flags to apply cuts at event level
932 +    vector<bool> eventPassedPreviousCuts;    //a vector to store cumulative cut descisions after each cut.
933      bool eventPassedAllCuts = true;
934 <
647 <    //apply trigger (true if none were specified)
648 <    eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
649 <
934 >    eventPassedAllCuts = eventPassedAllCuts && triggerDecision;    //apply trigger (true if none were specified)
935  
936      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
937  
938        //loop over all objects and count how many passed the cumulative selection up to this point
939        cut currentCut = currentChannel.cuts.at(currentCutIndex);
940        int numberPassing = 0;
941 +      int numberPassingPrev = 0;  // number of objects that pass cuts up to the previous one
942  
943        for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
944 <        if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
944 >        if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object).first) numberPassing++;
945        }
946 <      bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
946 >      bool cutDecision;
947 >      if (!currentCut.isVeto) {
948 >        cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
949 >      } else {
950 >        int prevCutIndex = currentCutIndex - 1;
951 >        if (prevCutIndex<0) {
952 >          numberPassingPrev = cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size();  // count all objects in collection if cut is the first
953 >        } else {
954 >          for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(prevCutIndex).size() ; object++){
955 >            if (cumulativeFlags.at(currentCut.inputCollection).at(prevCutIndex).at(object).first) {
956 >              numberPassingPrev++;
957 >              if (verbose_>1) clog << " object passed previous cut" << endl;  
958 >            }
959 >          }  
960 >        }  
961 >        int numberFailCut = numberPassingPrev - numberPassing;
962 >        cutDecision = evaluateComparison(numberFailCut,currentCut.eventComparativeOperator,currentCut.numberRequired);  
963 >      }
964 >
965        cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
966        eventPassedAllCuts = eventPassedAllCuts && cutDecision;
967 +      eventPassedPreviousCuts.push_back(eventPassedAllCuts);
968 +      if (verbose_>1) clog << " Event passed cuts up to cut index " << currentCutIndex << endl;  
969  
970      }
971 <
972 <    //     if(datasetType_ != "data") {
973 <    //       scaleFactor *= muonSFWeight_->at (chosenMuon ()->eta);
974 <    //       scaleFactor *= electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
975 <    //     }
976 <
971 >    //applying all appropriate scale factors
972 >    double scaleFactor = masterScaleFactor;
973 >    muonScaleFactor_ = electronScaleFactor_ = bTagScaleFactor_ = 1.0;
974 >
975 >    if(applyLeptonSF_ && datasetType_ != "data"){
976 >      //only apply SFs if we've cut on this object
977 >      if(find(objectsToCut.begin(),objectsToCut.end(),"muons") != objectsToCut.end ()){
978 >        flagPair muonFlags;
979 >        //get the last valid flags in the flag map
980 >        for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
981 >          if (cumulativeFlags.at("muons").at(i).size()){
982 >            muonFlags = cumulativeFlags.at("muons").at(i);
983 >            break;
984 >          }
985 >        }
986 >        //apply the weight for each of those objects
987 >        for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
988 >          if(!muonFlags.at(muonIndex).second) continue;
989 >          muonScaleFactor_ *= muonSFWeight_->at (muons->at(muonIndex).eta);
990 >        }
991 >      }
992 >      //only apply SFs if we've cut on this object
993 >      if(find(objectsToCut.begin(),objectsToCut.end(),"electrons") != objectsToCut.end ()){
994 >        flagPair electronFlags;
995 >        //get the last valid flags in the flag map
996 >        for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
997 >          if (cumulativeFlags.at("electrons").at(i).size()){
998 >            electronFlags = cumulativeFlags.at("electrons").at(i);
999 >            break;
1000 >          }
1001 >        }
1002 >        //apply the weight for each of those objects
1003 >        for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
1004 >          if(!electronFlags.at(electronIndex).second) continue;
1005 >          electronScaleFactor_ *= electronSFWeight_->at (electrons->at(electronIndex).eta, electrons->at(electronIndex).pt);
1006 >        }
1007 >      }
1008 >    }
1009 >    if(applyBtagSF_ && datasetType_ != "data"){
1010 >      //only apply SFs if we've cut on this object
1011 >      if(find(objectsToCut.begin(),objectsToCut.end(),"jets") != objectsToCut.end ()){
1012 >        flagPair jetFlags;
1013 >        vector<double> jetSFs;
1014 >        //get the last valid flags in the flag map
1015 >        for (int i = cumulativeFlags.at("jets").size() - 1; i >= 0; i--){
1016 >          if (cumulativeFlags.at("jets").at(i).size()){
1017 >            jetFlags = cumulativeFlags.at("jets").at(i);
1018 >            break;
1019 >          }
1020 >        }
1021 >        //apply the weight for each of those objects
1022 >        for (uint jetIndex = 0; jetIndex != jetFlags.size(); jetIndex++){
1023 >          if(!jetFlags.at(jetIndex).second) continue;
1024 >          double jetSFTmp = bTagSFWeight_->sflookup(jets->at(jetIndex).btagCombinedSecVertex, jets->at(jetIndex).pt, jets->at(jetIndex).flavour, jets->at(jetIndex).eta);
1025 >          jetSFs.push_back(jetSFTmp);
1026 >        }
1027 >        bTagScaleFactor_ *= bTagSFWeight_->weight( jetSFs, minBtag_);
1028 >      }
1029 >    }
1030 >    scaleFactor *= muonScaleFactor_;
1031 >    scaleFactor *= electronScaleFactor_;
1032 >    scaleFactor *= bTagScaleFactor_;
1033      cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
1034 +    if (verbose_>1) clog << " Scale factors applied:  "
1035 +                         << " muonScaleFactor_ = " << muonScaleFactor_
1036 +                         << ", electronScaleFactor_ = " << electronScaleFactor_
1037 +                         << ", bTagScaleFactor_ = " << bTagScaleFactor_
1038 +                         << ", total scale factor = " << scaleFactor
1039 +                         << endl;  
1040 +
1041 +
1042 +    if (printEventInfo_ && eventPassedAllCuts) {
1043 +      // Write information about event to screen, for testing purposes.
1044 +      clog << "Event passed all cuts in channel " <<  currentChannel.name
1045 +           << ":  run:lumi:evt=  "  << events->at(0).run
1046 +           << ":" << events->at(0).lumi
1047 +           << ":" << events->at(0).evt
1048 +           << endl;
1049 +      if (findEventsLog) {
1050 +        (*findEventsLog) << events->at(0).run
1051 +                         << ":" << events->at(0).lumi
1052 +                         << ":" << events->at(0).evt << endl;
1053 +      }
1054  
673
674
675    if(!eventPassedAllCuts)continue;
676
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;  
1055      }
1056  
1057  
1058      //filling histograms
1059 <    for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
1060 <      histogram currentHistogram = histograms.at(histogramIndex);
1059 >    for(uint currentCut = 0; currentCut != oneDHists_.at(currentChannelIndex).size(); currentCut++){//loop over all the cuts in each channel; if GetPlotsAfterEachCut_ is false, only the last cut will be used.
1060 >
1061 >      if (verbose_>2) clog << "  Filling histograms for currentcut = " << currentCut << endl;  
1062  
1063 <      if(currentHistogram.inputVariables.size() == 1){
1064 <        TH1D* histo;
1065 <        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
1066 <
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);
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);
1063 >      uint currentDir;
1064 >      if (!GetPlotsAfterEachCut_) { currentDir =  currentChannel.cuts.size() - oneDHists_.at(currentChannelIndex).size(); } //if GetPlotsAfterEachCut_ is false, set currentDir point to the last cut.
1065 >      else{
1066 >        currentDir = currentCut;
1067        }
787    }
1068  
1069  
1070 +      if(eventPassedPreviousCuts.at(currentDir)){
1071 +
1072 +        for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
1073 +          histogram currentHistogram = histograms.at(histogramIndex);
1074 +
1075 +          if (cumulativeFlags.count(currentHistogram.inputCollection) == 0) clog << "Error: no flags found for collection:  " << currentHistogram.inputCollection << ", will cause a seg fault" << endl;
1076 +
1077 +          if (verbose_>1) clog << " Filling histogram " << currentHistogram.name << " for collection " << currentHistogram.inputCollection << endl;  
1078 +
1079 +          if(currentHistogram.inputVariables.size() == 1){
1080 +            TH1D* histo;
1081 +            histo = oneDHists_.at(currentChannelIndex).at(currentCut).at(currentHistogram.name);
1082 +            if     (currentHistogram.inputCollection == "jets")            fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").at(currentDir),scaleFactor);
1083 +            else if(currentHistogram.inputCollection == "secondary jets")  fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("secondary jets").at(currentDir),scaleFactor);
1084 +            else if(currentHistogram.inputCollection == "secondary photons")  fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("secondary photons").at(currentDir),scaleFactor);
1085 +            else if(currentHistogram.inputCollection == "muons")           fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").at(currentDir),scaleFactor);
1086 +            else if(currentHistogram.inputCollection == "secondary muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").at(currentDir),scaleFactor);
1087 +            else if(currentHistogram.inputCollection == "secondary electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("secondary electrons").at(currentDir),scaleFactor);
1088 +            else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1089 +                                                                                           cumulativeFlags.at("muon-muon pairs").at(currentDir),scaleFactor);
1090 +            else if(currentHistogram.inputCollection == "muon-secondary muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1091 +                                                                                                     cumulativeFlags.at("muon-secondary muon pairs").at(currentDir),scaleFactor);
1092 +             else if(currentHistogram.inputCollection == "muon-secondary photon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1093 +                                                                                                     cumulativeFlags.at("muon-secondary photon pairs").at(currentDir),scaleFactor);
1094 +             else if(currentHistogram.inputCollection == "muon-secondary jet pairs") fill1DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1095 +                                                                                                     cumulativeFlags.at("muon-secondary jet pairs").at(currentDir),scaleFactor);
1096 +             else if(currentHistogram.inputCollection == "photon-secondary jet pairs") fill1DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1097 +                                                                                                     cumulativeFlags.at("photon-secondary jet pairs").at(currentDir),scaleFactor);
1098 +             else if(currentHistogram.inputCollection == "electron-secondary jet pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1099 +                                                                                                     cumulativeFlags.at("electron-secondary jet pairs").at(currentDir),scaleFactor);
1100 +
1101 +            else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").at(currentDir),scaleFactor);
1102 +            else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1103 +                                                                                                   cumulativeFlags.at("electron-electron pairs").at(currentDir),scaleFactor);
1104 +            else if(currentHistogram.inputCollection == "jet-jet pairs") fill1DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1105 +                                                                                         cumulativeFlags.at("jet-jet pairs").at(currentDir),scaleFactor);
1106 +             else if(currentHistogram.inputCollection == "jet-secondary jet pairs") fill1DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1107 +                                                                                                     cumulativeFlags.at("jet-secondary jet pairs").at(currentDir),scaleFactor);
1108 +
1109 +            else if(currentHistogram.inputCollection == "electron-secondary electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1110 +                                                                                                             cumulativeFlags.at("electron-secondary electron pairs").at(currentDir),scaleFactor);
1111 +            else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(),
1112 +                                                                                               cumulativeFlags.at("electron-muon pairs").at(currentDir),scaleFactor);
1113 +            else if(currentHistogram.inputCollection == "electron-jet pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),jets.product(),
1114 +                                                                                              cumulativeFlags.at("electron-jet pairs").at(currentDir),scaleFactor);
1115 +            else if(currentHistogram.inputCollection == "photon-jet pairs") fill1DHistogram(histo,currentHistogram, photons.product(),jets.product(),
1116 +                                                                                              cumulativeFlags.at("photon-jet pairs").at(currentDir),scaleFactor);
1117 +            else if(currentHistogram.inputCollection == "muon-jet pairs") fill1DHistogram(histo,currentHistogram, muons.product(),jets.product(),
1118 +                                                                                          cumulativeFlags.at("muon-jet pairs").at(currentDir),scaleFactor);
1119 +            else if(currentHistogram.inputCollection == "muon-event pairs") fill1DHistogram(histo,currentHistogram, muons.product(),events.product(),
1120 +                                                                                          cumulativeFlags.at("muon-event pairs").at(currentDir),scaleFactor);
1121 +            else if(currentHistogram.inputCollection == "met-jet pairs")  fill1DHistogram(histo,currentHistogram, mets.product(),jets.product(),
1122 +                                                                                          cumulativeFlags.at("met-jet pairs").at(currentDir),scaleFactor);
1123 +            else if(currentHistogram.inputCollection == "track-jet pairs")  fill1DHistogram(histo,currentHistogram,tracks.product(),jets.product(),
1124 +                                                                                          cumulativeFlags.at("track-jet pairs").at(currentDir),scaleFactor);
1125 +            else if(currentHistogram.inputCollection == "muon-photon pairs") fill1DHistogram(histo,currentHistogram, muons.product(),photons.product(),
1126 +                                                                                          cumulativeFlags.at("muon-photon pairs").at(currentDir),scaleFactor);
1127 +            else if(currentHistogram.inputCollection == "electron-photon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),photons.product(),
1128 +                                                                                          cumulativeFlags.at("electron-photon pairs").at(currentDir),scaleFactor);
1129 +            else if(currentHistogram.inputCollection == "electron-track pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),tracks.product(),
1130 +                                                                                                cumulativeFlags.at("electron-track pairs").at(currentDir),scaleFactor);
1131 +            else if(currentHistogram.inputCollection == "muon-track pairs") fill1DHistogram(histo,currentHistogram, muons.product(),tracks.product(),
1132 +                                                                                            cumulativeFlags.at("muon-track pairs").at(currentDir),scaleFactor);
1133 +            else if(currentHistogram.inputCollection == "muon-tau pairs") fill1DHistogram(histo,currentHistogram, muons.product(),taus.product(),
1134 +                                                                                          cumulativeFlags.at("muon-tau pairs").at(currentDir),scaleFactor);
1135 +            else if(currentHistogram.inputCollection == "tau-tau pairs") fill1DHistogram(histo,currentHistogram, taus.product(),taus.product(),
1136 +                                                                                         cumulativeFlags.at("tau-tau pairs").at(currentDir),scaleFactor);
1137 +            else if(currentHistogram.inputCollection == "tau-track pairs") fill1DHistogram(histo,currentHistogram, taus.product(),tracks.product(),
1138 +                                                                                           cumulativeFlags.at("tau-track pairs").at(currentDir),scaleFactor);
1139 +            else if(currentHistogram.inputCollection == "electron-trigobj pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),trigobjs.product(),
1140 +                                                                                                  cumulativeFlags.at("electron-trigobj pairs").at(currentDir),scaleFactor);
1141 +            else if(currentHistogram.inputCollection == "muon-trigobj pairs") fill1DHistogram(histo,currentHistogram, muons.product(),trigobjs.product(),
1142 +                                                                                              cumulativeFlags.at("muon-trigobj pairs").at(currentDir),scaleFactor);
1143 +
1144 +            else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").at(currentDir),scaleFactor);
1145 +            else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").at(currentDir),scaleFactor);
1146 +            else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").at(currentDir),scaleFactor);
1147 +            else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").at(currentDir),scaleFactor);
1148 +            else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").at(currentDir),scaleFactor);
1149 +            else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").at(currentDir),scaleFactor);
1150 +            else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").at(currentDir),scaleFactor);
1151 +            else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").at(currentDir),scaleFactor);
1152 +            else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").at(currentDir),scaleFactor);
1153 +            else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").at(currentDir),scaleFactor);
1154 +            else if(currentHistogram.inputCollection == "trigobjs") fill1DHistogram(histo,currentHistogram,trigobjs.product(),cumulativeFlags.at("trigobjs").at(currentDir),scaleFactor);
1155 +            else if(currentHistogram.inputCollection == "stops" && datasetType_ == "signalMC") fill1DHistogram(histo,currentHistogram,stops.product(),cumulativeFlags.at("stops").at(currentDir),scaleFactor);
1156 +          }
1157 +          else if(currentHistogram.inputVariables.size() == 2){
1158 +            TH2D* histo;
1159 +            histo = twoDHists_.at(currentChannelIndex).at(currentCut).at(currentHistogram.name);
1160 +            if     (currentHistogram.inputCollection == "jets")            fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").at(currentDir),scaleFactor);
1161 +            else if(currentHistogram.inputCollection == "secondary jets")  fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("secondary jets").at(currentDir),scaleFactor);
1162 +            else if(currentHistogram.inputCollection == "secondary photons")  fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("secondary photons").at(currentDir),scaleFactor);
1163 +            else if(currentHistogram.inputCollection == "muons")           fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").at(currentDir),scaleFactor);
1164 +            else if(currentHistogram.inputCollection == "secondary muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("secondary muons").at(currentDir),scaleFactor);
1165 +            else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1166 +                                                                                           cumulativeFlags.at("muon-muon pairs").at(currentDir),scaleFactor);
1167 +            else if(currentHistogram.inputCollection == "muon-secondary muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(),
1168 +                                                                                                     cumulativeFlags.at("muon-secondary muon pairs").at(currentDir),scaleFactor);
1169 +            else if(currentHistogram.inputCollection == "muon-secondary photon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1170 +                                                                                                     cumulativeFlags.at("muon-secondary photon pairs").at(currentDir),scaleFactor);
1171 +            else if(currentHistogram.inputCollection == "muon-secondary jet pairs") fill2DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1172 +                                                                                                     cumulativeFlags.at("muon-secondary jet pairs").at(currentDir),scaleFactor);
1173 +            else if(currentHistogram.inputCollection == "electron-secondary jet pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1174 +                                                                                                     cumulativeFlags.at("electron-secondary jet pairs").at(currentDir),scaleFactor);
1175 +            else if(currentHistogram.inputCollection == "photon-secondary jet pairs") fill2DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1176 +                                                                                                     cumulativeFlags.at("photon-secondary jet pairs").at(currentDir),scaleFactor);
1177 +            else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").at(currentDir),scaleFactor);
1178 +            else if(currentHistogram.inputCollection == "secondary electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("secondary electrons").at(currentDir),scaleFactor);
1179 +            else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1180 +                                                                                                   cumulativeFlags.at("electron-electron pairs").at(currentDir),scaleFactor);
1181 +            else if(currentHistogram.inputCollection == "jet-jet pairs") fill2DHistogram(histo,currentHistogram,jets.product(),jets.product(),
1182 +                                                                                         cumulativeFlags.at("jet-jet pairs").at(currentDir),scaleFactor);
1183 +            else if(currentHistogram.inputCollection == "electron-secondary electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),
1184 +                                                                                                             cumulativeFlags.at("electron-secondary electron pairs").at(currentDir),scaleFactor);
1185 +            else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(),
1186 +                                                                                               cumulativeFlags.at("electron-muon pairs").at(currentDir),scaleFactor);
1187 +            else if(currentHistogram.inputCollection == "electron-jet pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),jets.product(),
1188 +                                                                                              cumulativeFlags.at("electron-jet pairs").at(currentDir),scaleFactor);
1189 +            else if(currentHistogram.inputCollection == "electron-photon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),photons.product(),
1190 +                                                                                              cumulativeFlags.at("electron-photon pairs").at(currentDir),scaleFactor);
1191 +            else if(currentHistogram.inputCollection == "muon-jet pairs") fill2DHistogram(histo,currentHistogram,muons.product(),jets.product(),
1192 +                                                                                          cumulativeFlags.at("muon-jet pairs").at(currentDir),scaleFactor);
1193 +            else if(currentHistogram.inputCollection == "muon-event pairs") fill2DHistogram(histo,currentHistogram,muons.product(),events.product(),
1194 +                                                                                          cumulativeFlags.at("muon-event pairs").at(currentDir),scaleFactor);
1195 +            else if(currentHistogram.inputCollection == "met-jet pairs") fill2DHistogram(histo,currentHistogram,mets.product(),jets.product(),
1196 +                                                                                         cumulativeFlags.at("met-jet pairs").at(currentDir),scaleFactor);
1197 +            else if(currentHistogram.inputCollection == "track-jet pairs") fill2DHistogram(histo,currentHistogram,tracks.product(),jets.product(),
1198 +                                                                                         cumulativeFlags.at("track-jet pairs").at(currentDir),scaleFactor);
1199 +            else if(currentHistogram.inputCollection == "photon-jet pairs") fill2DHistogram(histo,currentHistogram,photons.product(),jets.product(),
1200 +                                                                                          cumulativeFlags.at("photon-jet pairs").at(currentDir),scaleFactor);
1201 +            else if(currentHistogram.inputCollection == "muon-photon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),photons.product(),
1202 +                                                                                          cumulativeFlags.at("muon-photon pairs").at(currentDir),scaleFactor);
1203 +            else if(currentHistogram.inputCollection == "electron-track pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),tracks.product(),
1204 +                                                                                                cumulativeFlags.at("electron-track pairs").at(currentDir),scaleFactor);
1205 +            else if(currentHistogram.inputCollection == "muon-track pairs") fill2DHistogram(histo,currentHistogram,muons.product(),tracks.product(),
1206 +                                                                                            cumulativeFlags.at("muon-track pairs").at(currentDir),scaleFactor);
1207 +            else if(currentHistogram.inputCollection == "muon-tau pairs") fill2DHistogram(histo,currentHistogram,muons.product(),taus.product(),
1208 +                                                                                          cumulativeFlags.at("muon-tau pairs").at(currentDir),scaleFactor);
1209 +            else if(currentHistogram.inputCollection == "tau-tau pairs") fill2DHistogram(histo,currentHistogram,taus.product(),taus.product(),
1210 +                                                                                         cumulativeFlags.at("tau-tau pairs").at(currentDir),scaleFactor);
1211 +            else if(currentHistogram.inputCollection == "tau-track pairs") fill2DHistogram(histo,currentHistogram,taus.product(),tracks.product(),
1212 +                                                                                           cumulativeFlags.at("tau-track pairs").at(currentDir),scaleFactor);
1213 +            else if(currentHistogram.inputCollection == "electron-trigobj pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),trigobjs.product(),
1214 +                                                                                                  cumulativeFlags.at("electron-trigobj pairs").at(currentDir),scaleFactor);
1215 +            else if(currentHistogram.inputCollection == "muon-trigobj pairs") fill2DHistogram(histo,currentHistogram,muons.product(),trigobjs.product(),
1216 +                                                                                              cumulativeFlags.at("muon-trigobj pairs").at(currentDir),scaleFactor);
1217 +            else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").at(currentDir),scaleFactor);
1218 +            else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").at(currentDir),scaleFactor);
1219 +            else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").at(currentDir),scaleFactor);
1220 +            else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").at(currentDir),scaleFactor);
1221 +            else if(currentHistogram.inputCollection == "track-event pairs") fill2DHistogram(histo,currentHistogram,tracks.product(),events.product(),
1222 +                                                                                             cumulativeFlags.at("track-event pairs").at(currentDir),scaleFactor);
1223 +            else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").at(currentDir),scaleFactor);
1224 +            else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").at(currentDir),scaleFactor);
1225 +            else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").at(currentDir),scaleFactor);
1226 +            else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").at(currentDir),scaleFactor);
1227 +            else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").at(currentDir),scaleFactor);
1228 +            else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").at(currentDir),scaleFactor);
1229 +            else if(currentHistogram.inputCollection == "trigobjs") fill2DHistogram(histo,currentHistogram,trigobjs.product(),cumulativeFlags.at("trigobjs").at(currentDir),scaleFactor);
1230 +            else if(currentHistogram.inputCollection == "stops" && datasetType_ == "signalMC") fill2DHistogram(histo,currentHistogram,stops.product(),cumulativeFlags.at("stops").at(currentDir),scaleFactor);
1231 +          }
1232  
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 = "";  
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";
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));  
814      string histoName = "num" + tempCurrentObject;
815
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();
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(),scaleFactor);
1233          }
827        else {
828          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,scaleFactor);
829        }
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      }
1234  
854    } // end for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
1235  
1236 <  } //end loop over channel
1236 >        //fills histograms with the sizes of collections
1237 >        for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
1238 >
1239 >          string currentObject = objectsToPlot.at(currentObjectIndex);
1240 >          string objectToPlot = "";
1241 >
1242 >          if (verbose_) clog << "Filling histogram of number of selected objects in collection: " << currentObject << endl;  
1243 >
1244 >          // Name of objectToPlot here must match the name specified in OSUAnalysis::OSUAnalysis().
1245 >          if(currentObject == "muon-muon pairs")                         objectToPlot = "dimuonPairs";
1246 >          else if(currentObject == "electron-electron pairs")            objectToPlot = "dielectronPairs";
1247 >          else if(currentObject == "electron-muon pairs")                objectToPlot = "electronMuonPairs";
1248 >          else if(currentObject == "electron-photon pairs")              objectToPlot = "electronPhotonPairs";
1249 >          else if(currentObject == "electron-jet pairs")                 objectToPlot = "electronJetPairs";
1250 >          else if(currentObject == "muon-jet pairs")                     objectToPlot = "muonJetPairs";
1251 >          else if(currentObject == "muon-event pairs")                   objectToPlot = "muonEventPairs";
1252 >          else if(currentObject == "muon-photon pairs")                  objectToPlot = "muonPhotonPairs";
1253 >          else if(currentObject == "photon-jet pairs")                   objectToPlot = "photonJetPairs";
1254 >          else if(currentObject == "met-jet pairs")                      objectToPlot = "metJetPairs";
1255 >          else if(currentObject == "track-jet pairs")                    objectToPlot = "trackJetPairs";
1256 >          else if(currentObject == "jet-jet pairs")                      objectToPlot = "dijetPairs";
1257 >          else if(currentObject == "jet-secondary jet pairs")            objectToPlot = "jetSecondaryJetPairs";
1258 >          else if(currentObject == "secondary jets")                     objectToPlot = "secondaryJets";
1259 >          else if(currentObject == "secondary photons")                  objectToPlot = "secondaryPhotons";
1260 >          else if(currentObject == "electron-track pairs")               objectToPlot = "electronTrackPairs";
1261 >          else if(currentObject == "muon-track pairs")                   objectToPlot = "muonTrackPairs";
1262 >          else if(currentObject == "muon-tau pairs")                     objectToPlot = "muonTauPairs";
1263 >          else if(currentObject == "tau-tau pairs")                      objectToPlot = "ditauPairs";
1264 >          else if(currentObject == "tau-track pairs")                    objectToPlot = "tauTrackPairs";
1265 >          else if(currentObject == "track-event pairs")                  objectToPlot = "trackEventPairs";
1266 >          else if(currentObject == "muon-secondary muon pairs")          objectToPlot = "muonSecondaryMuonPairs";
1267 >          else if(currentObject == "secondary muons")                    objectToPlot = "secondaryMuons";
1268 >          else if(currentObject == "muon-secondary jet pairs")           objectToPlot = "muonSecondaryJetPairs";
1269 >          else if(currentObject == "muon-secondary photon pairs")        objectToPlot = "muonSecondaryJetPairs";
1270 >          else if(currentObject == "electron-secondary jet pairs")       objectToPlot = "electronSecondaryJetPairs";
1271 >          else if(currentObject == "photon-secondary jet pairs")         objectToPlot = "photonSecondaryJetPairs";
1272 >          else if(currentObject == "electron-secondary electron pairs")  objectToPlot = "electronSecondaryElectronPairs";
1273 >          else if(currentObject == "secondary electrons")                objectToPlot = "secondaryElectrons";
1274 >          else if(currentObject == "electron-trigobj pairs")             objectToPlot = "electronTrigobjPairs";
1275 >          else if(currentObject == "muon-trigobj pairs")                 objectToPlot = "muonTrigobjPairs";
1276 >          else objectToPlot = currentObject;
1277 >
1278 >          string tempCurrentObject = objectToPlot;
1279 >          tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
1280 >          string histoName = "num" + tempCurrentObject;
1281 >
1282 >
1283 >          if(find(objectsToPlot.begin(), objectsToPlot.end(), currentObject) != objectsToPlot.end()) {
1284 >            flagPair lastCutFlags = cumulativeFlags.at(currentObject).at(currentDir);
1285 >            int numToPlot = 0;
1286 >            for (uint iObj = 0; iObj != lastCutFlags.size(); iObj++){  // loop over all the objects  
1287 >              if(lastCutFlags.at(iObj).second) {
1288 >                numToPlot++;  
1289 >                if (verbose_>3) clog << "   Found object " << iObj << " in collection " << currentObject << " to plot." << endl;
1290 >              }
1291 >            }
1292 >
1293 >            if(objectToPlot == "primaryvertexs"){
1294 >              oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
1295 >              oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
1296 >            }
1297 >            else {
1298 >              if (printEventInfo_) clog << "Number of selected " << objectToPlot << " to plot:  " << numToPlot << endl;  
1299 >              oneDHists_.at(currentChannelIndex).at(currentCut).at(histoName)->Fill(numToPlot,scaleFactor);
1300 >            }
1301 >          }
1302 >
1303 >        } // end for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++)
1304 >
1305  
1306 <  masterCutFlow_->fillCutFlow(scaleFactor);
1306 >      } // end if(eventPassedPreviousCuts.at(currentDir))
1307 >    } // end loop over cuts
1308  
1309 < } // end void OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
1309 >    if (!useEDMFormat_ && eventPassedAllCuts){
1310 >      // Assign BNTree variables
1311 >      for (uint iBranch = 0; iBranch != treeBranches_.size(); iBranch++) {
1312 >        BranchSpecs brSpecs = treeBranches_.at(iBranch);
1313 >        string coll = brSpecs.inputCollection;
1314 >        if (cumulativeFlags.count(coll) == 0) clog << "Error: no flags found for collection:  " << coll << ", will cause a seg fault" << endl;
1315 >
1316 >        if     (coll == "jets")                assignTreeBranch(brSpecs,jets.product(),          cumulativeFlags.at(coll).back());
1317 >        else if(coll == "secondary jets")      assignTreeBranch(brSpecs,jets.product(),          cumulativeFlags.at(coll).back());
1318 >        else if(coll == "secondary photons")      assignTreeBranch(brSpecs,photons.product(),          cumulativeFlags.at(coll).back());
1319 >        else if(coll == "muons")               assignTreeBranch(brSpecs,muons.product(),         cumulativeFlags.at(coll).back());
1320 >        else if(coll == "secondary muons")     assignTreeBranch(brSpecs,muons.product(),         cumulativeFlags.at(coll).back());
1321 >        else if(coll == "electrons")           assignTreeBranch(brSpecs,electrons.product(),     cumulativeFlags.at(coll).back());
1322 >        else if(coll == "secondary electrons") assignTreeBranch(brSpecs,electrons.product(),     cumulativeFlags.at(coll).back());
1323 >        else if(coll == "events")              assignTreeBranch(brSpecs,events.product(),        cumulativeFlags.at(coll).back());
1324 >        else if(coll == "taus")                assignTreeBranch(brSpecs,taus.product(),          cumulativeFlags.at(coll).back());
1325 >        else if(coll == "mets")                assignTreeBranch(brSpecs,mets.product(),          cumulativeFlags.at(coll).back());
1326 >        else if(coll == "tracks")              assignTreeBranch(brSpecs,tracks.product(),        cumulativeFlags.at(coll).back());
1327 >        else if(coll == "genjets")             assignTreeBranch(brSpecs,genjets.product(),       cumulativeFlags.at(coll).back());
1328 >        else if(coll == "mcparticles")         assignTreeBranch(brSpecs,mcparticles.product(),   cumulativeFlags.at(coll).back());
1329 >        else if(coll == "primaryvertexs")      assignTreeBranch(brSpecs,primaryvertexs.product(),cumulativeFlags.at(coll).back());
1330 >        else if(coll == "bxlumis")             assignTreeBranch(brSpecs,bxlumis.product(),       cumulativeFlags.at(coll).back());
1331 >        else if(coll == "photons")             assignTreeBranch(brSpecs,photons.product(),       cumulativeFlags.at(coll).back());
1332 >        else if(coll == "superclusters")       assignTreeBranch(brSpecs,superclusters.product(), cumulativeFlags.at(coll).back());
1333 >        else if(coll == "trigobjs")            assignTreeBranch(brSpecs,trigobjs.product(),      cumulativeFlags.at(coll).back());
1334 >        else if(coll == "stops"
1335 >                && datasetType_ == "signalMC") assignTreeBranch(brSpecs,stops.product(),         cumulativeFlags.at(coll).back());
1336 >      } // end loop over branches
1337 >
1338 >      if (!BNTrees_.at(currentChannelIndex)) { clog << "ERROR:  Undefined BNTree.  Will likely seg fault." << endl; }
1339 >      BNTrees_.at(currentChannelIndex)->Fill();  // only fill if event has passed cuts
1340 >    }
1341 >
1342 >    (*channelDecisions)[currentChannel.name] = eventPassedAllCuts;
1343 >
1344 >  } // end loop over channel
1345 >
1346 >  masterCutFlow_->fillCutFlow(masterScaleFactor);
1347 >
1348 >  event.put (channelDecisions, "channelDecisions");
1349 >
1350 >  isFirstEvent_ = false;  
1351 >
1352 >  if (verbose_) clog << "Finished OSUAnalysis::produce." << endl;  
1353 >
1354 > } // end void OSUAnalysis::produce (const edm::Event &event, const edm::EventSetup &setup)
1355  
1356  
1357  
# Line 872 | Line 1366 | OSUAnalysis::evaluateComparison (double
1366    else if(comparison == "==") return testValue == cutValue;
1367    else if(comparison == "=") return testValue == cutValue;
1368    else if(comparison == "!=") return testValue != cutValue;
1369 <  else {std::cout << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
1369 >  else {clog << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
1370  
1371   }
1372  
1373   bool
1374 < OSUAnalysis::evaluateTriggers (vector<string> triggersToTest, const BNtriggerCollection* triggerCollection){
1374 > OSUAnalysis::evaluateComparison (string testValue, string comparison, string cutValue){
1375 >
1376  
1377 +  if(comparison == ">")       return testValue >  cutValue;
1378 +  else if(comparison == ">=") return testValue >= cutValue;
1379 +  else if(comparison == "<")  return testValue <  cutValue;
1380 +  else if(comparison == "<=") return testValue <= cutValue;
1381 +  else if(comparison == "==") return testValue == cutValue;
1382 +  else if(comparison == "=")  return testValue == cutValue;
1383 +  else if(comparison == "!=") return testValue != cutValue;
1384 +  else {clog << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
1385 +
1386 + }
1387 +
1388 + bool
1389 + OSUAnalysis::evaluateTriggers (vector<string> triggersToTest, vector<string> triggersToVeto, const BNtriggerCollection* triggerCollection){
1390 +
1391 +  //initialize to false until a chosen trigger is passed
1392    bool triggerDecision = false;
1393  
1394 <  for (BNtriggerCollection::const_iterator trigger = triggerCollection->begin (); trigger != triggerCollection->end (); trigger++)
1395 <    {
1396 <      if(trigger->pass != 1) continue;
1397 <      for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++)
1398 <        {
1399 <          if(trigger->name.find(triggersToTest.at(triggerName))!=std::string::npos){
1400 <            triggerDecision = true;
1401 <          }
1402 <        }
1394 >  if (printAllTriggers_) clog << "Printing list of all available triggers (which this event may or may not pass):" << endl;
1395 >  //loop over all triggers defined in the event
1396 >  for (BNtriggerCollection::const_iterator trigger = triggerCollection->begin (); trigger != triggerCollection->end (); trigger++){
1397 >
1398 >    if (printAllTriggers_) clog << "   " << trigger->name << endl;
1399 >
1400 >    //we're only interested in triggers that actually passed
1401 >    if(trigger->pass != 1) continue;
1402 >
1403 >    //if the event passes one of the veto triggers, exit and return false
1404 >    for(uint triggerName = 0; triggerName != triggersToVeto.size(); triggerName++){
1405 >      if(trigger->name.find(triggersToVeto.at(triggerName))!=string::npos) return false;
1406      }
1407 <  return triggerDecision;
1407 >    //if the event passes one of the chosen triggers, set triggerDecision to true
1408 >    for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++){
1409 >      if(trigger->name.find(triggersToTest.at(triggerName))!=string::npos) triggerDecision = true;
1410 >    }
1411 >  }
1412 >
1413 >  printAllTriggers_ = false;  // only print triggers once, not every event
1414  
1415 +  //if none of the veto triggers fired:
1416 +  //return the OR of all the chosen triggers
1417 +  if (triggersToTest.size() != 0) return triggerDecision;
1418 +  //or if no triggers were defined return true
1419 +  else return true;
1420   }
1421  
1422 < std::vector<std::string>
1422 >
1423 > vector<string>
1424   OSUAnalysis::splitString (string inputString){
1425  
1426 <  std::stringstream stringStream(inputString);
1427 <  std::istream_iterator<std::string> begin(stringStream);
1428 <  std::istream_iterator<std::string> end;
1429 <  std::vector<std::string> stringVector(begin, end);
1426 >  stringstream stringStream(inputString);
1427 >  istream_iterator<string> begin(stringStream);
1428 >  istream_iterator<string> end;
1429 >  vector<string> stringVector(begin, end);
1430    return stringVector;
1431  
1432   }
1433  
1434  
1435   void OSUAnalysis::getTwoObjs(string tempInputCollection, string& obj1, string& obj2) {
1436 <  // Set two object strings from the tempInputCollection string,
1437 <  // For example, if tempInputCollection is "electron-muon pairs",
1438 <  // then obj1 = "electrons" and obj2 = "muons".  
1439 <  // Note that the objects have an "s" appended.  
1436 >  // Set two object strings from the tempInputCollection string,
1437 >  // For example, if tempInputCollection is "electron-muon pairs",
1438 >  // then obj1 = "electrons" and obj2 = "muons".
1439 >  // Note that the objects have an "s" appended.
1440  
1441    int dashIndex = tempInputCollection.find("-");
1442    int spaceIndex = tempInputCollection.find_last_of(" ");
1443    int secondWordLength = spaceIndex - dashIndex;
1444 <  obj1 = tempInputCollection.substr(0,dashIndex) + "s";  
1444 >  obj1 = tempInputCollection.substr(0,dashIndex) + "s";
1445    obj2 = tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s";
1446 <  
1447 < }
1446 >
1447 > }
1448  
1449  
1450   string OSUAnalysis::getObjToGet(string obj) {
1451    // Return the string corresponding to the object to get for the given obj string.
1452 <  // Right now this only handles the case in which obj contains "secondary",
1453 <  // e.g, "secondary muons".  
1454 <  // 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);
1452 >  // Right now this only handles the case in which obj contains "secondary",
1453 >  // e.g, "secondary muons".
1454 >  // Note that "s" is NOT appended.
1455  
1456 < }  
1456 >  if (obj.find("secondary")==string::npos) return obj;  // "secondary" is not found
1457 >  int firstSpaceIndex = obj.find_first_of(" ");
1458 >  return obj.substr(firstSpaceIndex+1,obj.length()-1);
1459  
1460 + }
1461  
1462  
1463 + //!jet valueLookup
1464   double
1465 < OSUAnalysis::valueLookup (const BNjet* object, string variable, string function){
1465 > OSUAnalysis::valueLookup (const BNjet* object, string variable, string function, string &stringValue){
1466  
1467    double value = 0.0;
1468    if(variable == "energy") value = object->energy;
# Line 1056 | Line 1581 | OSUAnalysis::valueLookup (const BNjet* o
1581    else if(variable == "puJetId_loose_simple") value = object->puJetId_loose_simple;
1582    else if(variable == "puJetId_loose_cutbased") value = object->puJetId_loose_cutbased;
1583  
1584 +  //user defined variable
1585 +  else if(variable == "disappTrkLeadingJetID") {
1586 +    value = object->pt > 110
1587 +      && fabs(object->eta) < 2.4
1588 +      && object->chargedHadronEnergyFraction > 0.2
1589 +      && object->neutralHadronEnergyFraction < 0.7
1590 +      && object->chargedEmEnergyFraction < 0.5
1591 +      && object->neutralEmEnergyFraction < 0.7;
1592 +  }
1593  
1594 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1594 >  else if(variable == "disappTrkSubLeadingJetID") {
1595 >    value = object->pt > 30
1596 >      && fabs(object->eta) < 4.5
1597 >      && object->neutralHadronEnergyFraction < 0.7
1598 >      && object->chargedEmEnergyFraction < 0.5;
1599 >  }
1600 >
1601 >  else if(variable == "dPhiMet") {
1602 >    if (const BNmet *met = chosenMET ()) {
1603 >      value = deltaPhi (object->phi, met->phi);
1604 >    } else value = -999;
1605 >  }
1606 >
1607 >
1608 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1609  
1610    value = applyFunction(function, value);
1611  
1612    return value;
1613 < }
1066 <
1613 > } // end jet valueLookup
1614  
1615  
1616 + //!muon valueLookup
1617   double
1618 < OSUAnalysis::valueLookup (const BNmuon* object, string variable, string function){
1618 > OSUAnalysis::valueLookup (const BNmuon* object, string variable, string function, string &stringValue){
1619  
1620    double value = 0.0;
1621    if(variable == "energy") value = object->energy;
# Line 1159 | Line 1707 | OSUAnalysis::valueLookup (const BNmuon*
1707    else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1708    else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1709    else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1710 +  else if(variable == "chargedHadronIsoDR04") value = object->chargedHadronIsoDR04;
1711 +  else if(variable == "neutralHadronIsoDR04") value = object->neutralHadronIsoDR04;
1712 +  else if(variable == "photonIsoDR04") value = object->photonIsoDR04;
1713 +  else if(variable == "puChargedHadronIsoDR04") value = object->puChargedHadronIsoDR04;
1714    else if(variable == "rhoPrime") value = object->rhoPrime;
1715    else if(variable == "AEffDr03") value = object->AEffDr03;
1716    else if(variable == "AEffDr04") value = object->AEffDr04;
# Line 1166 | Line 1718 | OSUAnalysis::valueLookup (const BNmuon*
1718    else if(variable == "pfIsoR03SumNeutralHadronEt") value = object->pfIsoR03SumNeutralHadronEt;
1719    else if(variable == "pfIsoR03SumPhotonEt") value = object->pfIsoR03SumPhotonEt;
1720    else if(variable == "pfIsoR03SumPUPt") value = object->pfIsoR03SumPUPt;
1721 +  else if(variable == "relpfIsoR04SumExceptChargedHad") value = (object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)/ object->pt;
1722 +  else if(variable == "relpfIsoR04SumChargedHadronPt") value = object->pfIsoR04SumChargedHadronPt/object->pt;
1723 +  else if(variable == "relpfIsoR04SumNeutralHadronEt") value = object->pfIsoR04SumNeutralHadronEt/object->pt;
1724 +  else if(variable == "relpfIsoR04SumPhotonEt") value = object->pfIsoR04SumPhotonEt/object->pt;
1725 +  else if(variable == "relpfIsoR04SumPUPt") value = object->pfIsoR04SumPUPt/object->pt;
1726    else if(variable == "pfIsoR04SumChargedHadronPt") value = object->pfIsoR04SumChargedHadronPt;
1727    else if(variable == "pfIsoR04SumNeutralHadronEt") value = object->pfIsoR04SumNeutralHadronEt;
1728    else if(variable == "pfIsoR04SumPhotonEt") value = object->pfIsoR04SumPhotonEt;
# Line 1222 | Line 1779 | OSUAnalysis::valueLookup (const BNmuon*
1779    else if(variable == "time_ndof") value = object->time_ndof;
1780  
1781    //user-defined variables
1782 +  else if(variable == "looseID") {
1783 +    value = object->pt > 10 &&
1784 +      (object->isGlobalMuon  > 0 ||
1785 +       object->isTrackerMuon > 0);
1786 +  }
1787    else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1788    else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1789    else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1790    else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1791 +  else if(variable == "relPFdBetaIsoPseudo") value = (object->pfIsoR04SumChargedHadronPt + object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt) / object->pt;
1792    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1793    else if(variable == "metMT") {
1794 <    const BNmet *met = chosenMET ();
1232 <    if (met)
1794 >    if (const BNmet *met = chosenMET ())
1795        {
1796 <        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1797 <          p2 (met->px, met->py, 0.0, met->pt);
1236 <
1237 <        value = (p1 + p2).Mt ();
1796 >        double dPhi = deltaPhi (object->phi, met->phi);
1797 >        value = sqrt (2 * object->pt * met->pt * (1 - cos (dPhi)));
1798        }
1799      else
1800        value = -999;
# Line 1352 | Line 1912 | OSUAnalysis::valueLookup (const BNmuon*
1912      value = object->isGlobalMuon > 0                \
1913        && object->isPFMuon > 0                        \
1914        && object->normalizedChi2 < 10                \
1915 <                                  && object->numberOfValidMuonHits > 0        \
1915 >      && object->numberOfValidMuonHits > 0        \
1916        && object->numberOfMatchedStations > 1        \
1917        && fabs(object->correctedD0Vertex) < 0.2        \
1918        && fabs(object->correctedDZ) < 0.5        \
# Line 1361 | Line 1921 | OSUAnalysis::valueLookup (const BNmuon*
1921    }
1922    else if(variable == "tightIDdisplaced"){
1923      value = object->isGlobalMuon > 0                \
1924 +      && object->isPFMuon > 0                        \
1925        && object->normalizedChi2 < 10                \
1926 <                                  && object->numberOfValidMuonHits > 0        \
1926 >      && object->numberOfValidMuonHits > 0        \
1927        && object->numberOfMatchedStations > 1        \
1928        && object->numberOfValidPixelHits > 0        \
1929        && object->numberOfLayersWithMeasurement > 5;
1930    }
1931  
1932 <  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1932 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1933  
1934    else if(variable == "genMatchedPdgId"){
1935      int index = getGenMatchedParticleIndex(object);
# Line 1388 | Line 1949 | OSUAnalysis::valueLookup (const BNmuon*
1949    }
1950    else if(variable == "genMatchedMotherIdReverse"){
1951      int index = getGenMatchedParticleIndex(object);
1952 <    if(index == -1) value = 23;
1953 <    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
1952 >    if(index == -1) value = 24;
1953 >    else value = 24 - getPdgIdBinValue(mcparticles->at(index).motherId);
1954    }
1955    else if(variable == "genMatchedGrandmotherId"){
1956      int index = getGenMatchedParticleIndex(object);
# Line 1401 | Line 1962 | OSUAnalysis::valueLookup (const BNmuon*
1962      }
1963      else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1964    }
1965 +  else if(variable == "genMatchedGrandmotherIdReverse"){
1966 +    int index = getGenMatchedParticleIndex(object);
1967 +    if(index == -1) value = 24;
1968 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1969 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1970 +      if(motherIndex == -1) value = 24;
1971 +      else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1972 +    }
1973 +    else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1974 +  }
1975 +  else if(variable == "pfMuonsFromVertex"){
1976 +    double d0Error, dzError;
1977 +
1978 +    d0Error = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1979 +    dzError = hypot (object->tkDZerr, chosenVertex ()->zError);
1980 +    value = fabs (object->correctedD0Vertex) > 0.2 || fabs (object->correctedDZ) > 0.5
1981 +      || fabs (object->correctedD0Vertex / d0Error) > 99.0
1982 +      || fabs (object->correctedDZ / dzError) > 99.0;
1983 +    value = (object->isStandAloneMuon && !object->isTrackerMuon && !object->isGlobalMuon) || !value;
1984 +  }
1985  
1986  
1987  
1988 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1988 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1989  
1990    value = applyFunction(function, value);
1991  
1992    return value;
1993 < }
1993 > } // end muon valueLookup
1994  
1995  
1996 + //!electron valueLookup
1997   double
1998 < OSUAnalysis::valueLookup (const BNelectron* object, string variable, string function){
1998 > OSUAnalysis::valueLookup (const BNelectron* object, string variable, string function, string &stringValue){
1999  
2000    double value = 0.0;
2001    if(variable == "energy") value = object->energy;
# Line 1578 | Line 2160 | OSUAnalysis::valueLookup (const BNelectr
2160    else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
2161    else if(variable == "detIso") value = (object->trackIso) / object->pt;
2162    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
2163 +  else if(variable == "relPFrhoIsoEB") value = object->isEB ? ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt : -999;
2164 +  else if(variable == "relPFrhoIsoEE") value = object->isEE ? ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt : -999;
2165    else if(variable == "metMT") {
2166 <    const BNmet *met = chosenMET ();
1583 <    if (met)
2166 >    if (const BNmet *met = chosenMET ())
2167        {
2168 <        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
2169 <          p2 (met->px, met->py, 0.0, met->pt);
1587 <
1588 <        value = (p1 + p2).Mt ();
2168 >        double dPhi = deltaPhi (object->phi, met->phi);
2169 >        value = sqrt (2 * object->pt * met->pt * (1 - cos (dPhi)));
2170        }
2171      else
2172        value = -999;
# Line 1652 | Line 2233 | OSUAnalysis::valueLookup (const BNelectr
2233      else value = -999;
2234    }
2235  
2236 +  else if(variable == "looseID"){
2237 +    if (object->isEB)
2238 +      {
2239 +        value = fabs(object->delEtaIn) < 0.007 \
2240 +          && fabs (object->delPhiIn) < 0.15 \
2241 +          && object->scSigmaIEtaIEta < 0.01 \
2242 +          && object->hadOverEm < 0.12 \
2243 +          && fabs (object->correctedD0Vertex) < 0.02 \
2244 +          && fabs (object->correctedDZ) < 0.2 \
2245 +          && object->absInvEMinusInvPin < 0.05 \
2246 +          && object->passConvVeto;
2247 +      }
2248 +    else
2249 +      {
2250 +        value = fabs(object->delEtaIn) < 0.009 \
2251 +          && fabs (object->delPhiIn) < 0.10 \
2252 +          && object->scSigmaIEtaIEta < 0.03 \
2253 +          && object->hadOverEm < 0.10 \
2254 +          && fabs (object->correctedD0Vertex) < 0.02 \
2255 +          && fabs (object->correctedDZ) < 0.2 \
2256 +          && object->absInvEMinusInvPin < 0.05 \
2257 +          && object->passConvVeto;
2258 +      }
2259 +  }
2260 +
2261    else if(variable == "tightID"){
2262      if (object->isEB)
2263        {
# Line 1677 | Line 2283 | OSUAnalysis::valueLookup (const BNelectr
2283        }
2284    }
2285  
2286 +  else if(variable == "looseID_MVA"){
2287 +    value = object->pt > 10
2288 +      && object->mvaNonTrigV0 > 0;
2289 +      }
2290    else if(variable == "correctedD0VertexInEBPositiveCharge"){
2291      if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
2292      else value = -999;
# Line 1750 | Line 2360 | OSUAnalysis::valueLookup (const BNelectr
2360    }
2361  
2362  
2363 <  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2363 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2364  
2365    else if(variable == "genMatchedPdgId"){
2366      int index = getGenMatchedParticleIndex(object);
# Line 1771 | Line 2381 | OSUAnalysis::valueLookup (const BNelectr
2381    }
2382    else if(variable == "genMatchedMotherIdReverse"){
2383      int index = getGenMatchedParticleIndex(object);
2384 <    if(index == -1) value = 23;
2385 <    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2384 >    if(index == -1) value = 24;
2385 >    else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2386    }
2387    else if(variable == "genMatchedGrandmotherId"){
2388      int index = getGenMatchedParticleIndex(object);
# Line 1784 | Line 2394 | OSUAnalysis::valueLookup (const BNelectr
2394      }
2395      else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2396    }
2397 +  else if(variable == "genMatchedGrandmotherIdReverse"){
2398 +    int index = getGenMatchedParticleIndex(object);
2399 +    if(index == -1) value = 24;
2400 +    else if(fabs(mcparticles->at(index).motherId) == 15){
2401 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2402 +      if(motherIndex == -1) value = 24;
2403 +      else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2404 +    }
2405 +    else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2406 +  }
2407 +  else if(variable == "pfElectronsFromVertex"){
2408 +    double d0Error, dzError;
2409 +
2410 +    d0Error = hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
2411 +    dzError = hypot (object->tkDZerr, chosenVertex ()->zError);
2412 +    value = fabs (object->correctedD0Vertex) > 0.2 || fabs (object->correctedDZ) > 0.5
2413 +      || fabs (object->correctedD0Vertex / d0Error) > 99.0
2414 +      || fabs (object->correctedDZ / dzError) > 99.0;
2415 +    value = !value;
2416 +  }
2417  
2418  
2419  
2420 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2420 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2421  
2422    value = applyFunction(function, value);
2423  
2424    return value;
2425 < }
2425 > } // end electron valueLookup
2426  
2427  
2428 + //!event valueLookup
2429   double
2430 < OSUAnalysis::valueLookup (const BNevent* object, string variable, string function){
2430 > OSUAnalysis::valueLookup (const BNevent* object, string variable, string function, string &stringValue){
2431  
2432    double value = 0.0;
2433  
# Line 1804 | Line 2435 | OSUAnalysis::valueLookup (const BNevent*
2435    else if(variable == "pthat") value = object->pthat;
2436    else if(variable == "qScale") value = object->qScale;
2437    else if(variable == "alphaQCD") value = object->alphaQCD;
2438 +  else if(variable == "Ht") value = getHt(jets.product());
2439    else if(variable == "alphaQED") value = object->alphaQED;
2440    else if(variable == "scalePDF") value = object->scalePDF;
2441    else if(variable == "x1") value = object->x1;
# Line 1865 | Line 2497 | OSUAnalysis::valueLookup (const BNevent*
2497    else if(variable == "id2") value = object->id2;
2498    else if(variable == "evt") value = object->evt;
2499    else if(variable == "puScaleFactor"){
2500 <    if(datasetType_ != "data")
2500 >    if(doPileupReweighting_ && datasetType_ != "data")
2501        value = puWeight_->at (events->at (0).numTruePV);
2502      else
2503        value = 1.0;
2504    }
2505 <  else if(variable == "muonScaleFactor"){
2506 <    if(datasetType_ != "data")
2507 <      //      value = muonSFWeight_->at (chosenMuon ()->eta);
2508 <      value = 1.0;
2509 <    else
2510 <      value = 1.0;
2511 <  }
2512 <  else if(variable == "electronScaleFactor"){
2513 <    if(datasetType_ != "data")
2514 <      //      value = electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
2515 <      value = 1.0;
2516 <    else
2517 <      value = 1.0;
2505 >  else if(variable == "muonScaleFactor") value = muonScaleFactor_;
2506 >  else if(variable == "electronScaleFactor") value = electronScaleFactor_;
2507 >  else if(variable == "stopCTauScaleFactor") value = stopCTauScaleFactor_;
2508 >  else if(variable == "bTagScaleFactor") value = bTagScaleFactor_;
2509 >  else if(variable == "ht") value = chosenHT ();
2510 >  else if(variable == "leadMuPairInvMass"){
2511 >    pair<const BNmuon *, const BNmuon *> muPair = leadMuonPair ();
2512 >    TLorentzVector p0 (muPair.first->px, muPair.first->py, muPair.first->pz, muPair.first->energy),
2513 >                   p1 (muPair.second->px, muPair.second->py, muPair.second->pz, muPair.second->energy);
2514 >    value = (p0 + p1).M ();
2515 >  }
2516 >  else if(variable == "leadMuPairPt"){
2517 >    pair<const BNmuon *, const BNmuon *> muPair = leadMuonPair ();
2518 >    TVector2 pt0 (muPair.first->px, muPair.first->py),
2519 >             pt1 (muPair.second->px, muPair.second->py);
2520 >    pt0 += pt1;
2521 >    value = pt0.Mod ();
2522 >  }
2523 >  else if(variable == "leadElPairInvMass"){
2524 >    pair<const BNelectron *, const BNelectron *> muPair = leadElectronPair ();
2525 >    TLorentzVector p0 (muPair.first->px, muPair.first->py, muPair.first->pz, muPair.first->energy),
2526 >                   p1 (muPair.second->px, muPair.second->py, muPair.second->pz, muPair.second->energy);
2527 >    value = (p0 + p1).M ();
2528 >  }
2529 >  else if(variable == "leadElPairPt"){
2530 >    pair<const BNelectron *, const BNelectron *> muPair = leadElectronPair ();
2531 >    TVector2 pt0 (muPair.first->px, muPair.first->py),
2532 >             pt1 (muPair.second->px, muPair.second->py);
2533 >    pt0 += pt1;
2534 >    value = pt0.Mod ();
2535    }
2536 <  else if(variable == "cTauScaleFactor")
1888 <    value = cTauScaleFactor_;
1889 <
1890 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2536 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2537  
2538    value = applyFunction(function, value);
2539  
2540    return value;
2541 < }
2541 > } // end event valueLookup
2542  
2543 +
2544 + //!tau valueLookup
2545   double
2546 < OSUAnalysis::valueLookup (const BNtau* object, string variable, string function){
2546 > OSUAnalysis::valueLookup (const BNtau* object, string variable, string function, string &stringValue){
2547  
2548    double value = 0.0;
2549  
# Line 1939 | Line 2587 | OSUAnalysis::valueLookup (const BNtau* o
2587    else if(variable == "HPSdecayModeFinding") value = object->HPSdecayModeFinding;
2588    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
2589  
2590 +  else if (variable == "looseHadronicID") {
2591 +    value = object->pt > 10
2592 +      && object->eta < 2.3
2593 +      && object->HPSbyLooseCombinedIsolationDeltaBetaCorr > 0
2594 +      && object->HPSdecayModeFinding > 0
2595 +      && object->HPSagainstElectronLoose > 0
2596 +      && object->HPSagainstMuonTight > 0;
2597 +  }
2598  
2599 <
1944 <  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2599 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2600  
2601    else if(variable == "genMatchedPdgId"){
2602      int index = getGenMatchedParticleIndex(object);
# Line 1961 | Line 2616 | OSUAnalysis::valueLookup (const BNtau* o
2616    }
2617    else if(variable == "genMatchedMotherIdReverse"){
2618      int index = getGenMatchedParticleIndex(object);
2619 <    if(index == -1) value = 23;
2620 <    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2619 >    if(index == -1) value = 24;
2620 >    else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2621    }
2622    else if(variable == "genMatchedGrandmotherId"){
2623      int index = getGenMatchedParticleIndex(object);
# Line 1974 | Line 2629 | OSUAnalysis::valueLookup (const BNtau* o
2629      }
2630      else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2631    }
2632 +  else if(variable == "genMatchedGrandmotherIdReverse"){
2633 +    int index = getGenMatchedParticleIndex(object);
2634 +    if(index == -1) value = 24;
2635 +    else if(fabs(mcparticles->at(index).motherId) == 15){
2636 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2637 +      if(motherIndex == -1) value = 24;
2638 +      else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2639 +    }
2640 +    else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2641 +  }
2642  
2643  
2644 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2644 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2645  
2646    value = applyFunction(function, value);
2647  
2648    return value;
2649 < }
2649 > } // end tau valueLookup
2650 >
2651  
2652 + //!met valueLookup
2653   double
2654 < OSUAnalysis::valueLookup (const BNmet* object, string variable, string function){
2654 > OSUAnalysis::valueLookup (const BNmet* object, string variable, string function, string &stringValue){
2655  
2656    double value = 0.0;
2657  
# Line 2048 | Line 2715 | OSUAnalysis::valueLookup (const BNmet* o
2715    else if(variable == "pfT1jet10pt") value = object->pfT1jet10pt;
2716    else if(variable == "pfT1jet10phi") value = object->pfT1jet10phi;
2717  
2718 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2718 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2719  
2720    value = applyFunction(function, value);
2721  
2722    return value;
2723 < }
2723 > } // end met valueLookup
2724 >
2725  
2726 + //!track valueLookup
2727   double
2728 < OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){
2728 > OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function, string &stringValue){
2729  
2730    double value = 0.0;
2731    double pMag = sqrt(object->pt * object->pt +
2732 <                     object->pz * object->pz);
2733 <  
2732 >                     object->pz * object->pz);
2733 >
2734    if(variable == "pt") value = object->pt;
2735    else if(variable == "px") value = object->px;
2736    else if(variable == "py") value = object->py;
# Line 2080 | Line 2749 | OSUAnalysis::valueLookup (const BNtrack*
2749    else if(variable == "numValidHits") value = object->numValidHits;
2750    else if(variable == "isHighPurity") value = object->isHighPurity;
2751  
2083
2752    //additional BNs info for disappTrks
2753 <  else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
2754 <  else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
2755 <  else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
2756 <  else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
2757 <  else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
2758 <  else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
2759 <  else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
2760 <  else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
2093 <  else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
2753 >  else if(variable == "caloEMDeltaRp3")     value = object->caloEMDeltaRp3;
2754 >  else if(variable == "caloHadDeltaRp3")    value = object->caloHadDeltaRp3;
2755 >  else if(variable == "caloEMDeltaRp4")     value = object->caloEMDeltaRp4;
2756 >  else if(variable == "caloHadDeltaRp4")    value = object->caloHadDeltaRp4;
2757 >  else if(variable == "caloEMDeltaRp5")     value = object->caloEMDeltaRp5;
2758 >  else if(variable == "caloHadDeltaRp5")    value = object->caloHadDeltaRp5;
2759 >  else if(variable == "nHitsMissingOuter")  value = object->nHitsMissingOuter;
2760 >  else if(variable == "nHitsMissingInner")  value = object->nHitsMissingInner;
2761    else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
2762 <  
2762 >  else if(variable == "depTrkRp3")          value = object->depTrkRp3;
2763 >  else if(variable == "trkRelIsoRp3")       value = (object->depTrkRp3 - object->pt) / object->pt;
2764 >  else if(variable == "trkRelIsoRp5")       value = (object->depTrkRp5 - object->pt) / object->pt;
2765 >  else if(variable == "depEcalRp3")         value = object->depEcalRp3;
2766 >  else if(variable == "depHcalRp3")         value = object->depHcalRp3;
2767 >  else if(variable == "depHoRp3")           value = object->depHoRp3;
2768 >  else if(variable == "nTracksRp3")         value = object->nTracksRp3;
2769 >  else if(variable == "trackerVetoPtRp3")   value = object->trackerVetoPtRp3;
2770 >  else if(variable == "emVetoEtRp3")        value = object->emVetoEtRp3;
2771 >  else if(variable == "hadVetoEtRp3")       value = object->hadVetoEtRp3;
2772 >  else if(variable == "hoVetoEtRp3")        value = object->hoVetoEtRp3;
2773 >  else if(variable == "depTrkRp5")          value = object->depTrkRp5;
2774 >  else if(variable == "depEcalRp5")         value = object->depEcalRp5;
2775 >  else if(variable == "depHcalRp5")         value = object->depHcalRp5;
2776 >  else if(variable == "depHoRp5")           value = object->depHoRp5;
2777 >  else if(variable == "nTracksRp5")         value = object->nTracksRp5;
2778 >  else if(variable == "trackerVetoPtRp5")   value = object->trackerVetoPtRp5;
2779 >  else if(variable == "emVetoEtRp5")        value = object->emVetoEtRp5;
2780 >  else if(variable == "hadVetoEtRp5")       value = object->hadVetoEtRp5;
2781 >  else if(variable == "hoVetoEtRp5")        value = object->hoVetoEtRp5;
2782  
2783    //user defined variables
2784    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;
2785    else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
2786 +  else if(variable == "depTrkRp5MinusPt"){
2787 +    if ( (object->depTrkRp5 - object->pt) < 0 ) {
2788 + //       clog << "Warning:  found track with depTrkRp5 < pt:  depTrkRp5=" << object->depTrkRp5
2789 + //         << "; pt=" << object->pt
2790 + //         << "; object->depTrkRp5 - object->pt = " << object->depTrkRp5 - object->pt
2791 + //         << endl;  
2792 +           value = 0;
2793 +         }
2794 +         else value =  (object->depTrkRp5 - object->pt);
2795 +  }
2796 +  else if(variable == "depTrkRp3MinusPt"){
2797 +    if ( (object->depTrkRp3 - object->pt) < 0 ) {
2798 +      value = 0;
2799 +    }
2800 +    else value =  (object->depTrkRp3 - object->pt);
2801 +  }
2802 +
2803 +  else if(variable == "dPhiMet") {
2804 +    if (const BNmet *met = chosenMET ()) {
2805 +      value = deltaPhi (object->phi, met->phi);
2806 +    } else value = -999;
2807 +  }
2808 +  
2809 +  
2810    else if(variable == "caloTotDeltaRp5")            value =  (object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
2811    else if(variable == "caloTotDeltaRp5ByP")         value = ((object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
2812 <  else if(variable == "caloTotDeltaRp5RhoCorr")     value = getTrkCaloTotRhoCorr(object);  
2813 <  else if(variable == "caloTotDeltaRp5ByPRhoCorr")  value = getTrkCaloTotRhoCorr(object) / pMag;  
2812 >  else if(variable == "caloTotDeltaRp5RhoCorr")     value = getTrkCaloTotRhoCorr(object);
2813 >  else if(variable == "caloTotDeltaRp5ByPRhoCorr")  value = getTrkCaloTotRhoCorr(object) / pMag;
2814 >  else if(variable == "depTrkRp5RhoCorr")           value = getTrkDepTrkRp5RhoCorr(object);
2815 >  else if(variable == "depTrkRp3RhoCorr")           value = getTrkDepTrkRp3RhoCorr(object);
2816 >
2817 >  else if(variable == "depTrkRp5MinusPtRhoCorr")    {
2818 >    if ( (getTrkDepTrkRp5RhoCorr(object) - object->pt) < 0 ) value = 0;
2819 >    else {value = (getTrkDepTrkRp5RhoCorr(object) - object->pt );}
2820 >  }
2821 >  
2822 >  else if(variable == "depTrkRp3MinusPtRhoCorr")    
2823 >    {
2824 >      if ( (getTrkDepTrkRp3RhoCorr(object) - object->pt) < 0 ) value = 0;
2825 >      else {value = (getTrkDepTrkRp3RhoCorr(object) - object->pt );}
2826 >    }
2827 >
2828    else if(variable == "isIso")                      value = getTrkIsIso(object, tracks.product());
2829    else if(variable == "isMatchedDeadEcal")          value = getTrkIsMatchedDeadEcal(object);
2830    else if(variable == "ptErrorByPt")                value = (object->ptError/object->pt);
2831    else if(variable == "ptError")                    value = object->ptError;
2832    else if(variable == "ptRes")                      value = getTrkPtRes(object);
2833 <  else if (variable == "d0wrtPV"){      
2834 <    double vx = object->vx - chosenVertex ()->x,        
2835 <      vy = object->vy - chosenVertex ()->y,      
2836 <      px = object->px,  
2837 <      py = object->py,  
2838 <      pt = object->pt;  
2839 <    value = (-vx * py + vy * px) / pt;  
2840 <  }      
2841 <  else if (variable == "dZwrtPV"){      
2842 <    double vx = object->vx - chosenVertex ()->x,        
2843 <      vy = object->vy - chosenVertex ()->y,      
2844 <      vz = object->vz - chosenVertex ()->z,      
2845 <      px = object->px,  
2846 <      py = object->py,  
2847 <      pz = object->pz,  
2848 <      pt = object->pt;  
2849 <    value = vz - (vx * px + vy * py)/pt * (pz/pt);      
2850 <  }    
2851 <  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2833 >  else if (variable == "d0wrtPV"){
2834 >    double vx = object->vx - chosenVertex ()->x,
2835 >      vy = object->vy - chosenVertex ()->y,
2836 >      px = object->px,
2837 >      py = object->py,
2838 >      pt = object->pt;
2839 >    value = (-vx * py + vy * px) / pt;
2840 >  }
2841 >  else if (variable == "dZwrtPV"){
2842 >    double vx = object->vx - chosenVertex ()->x,
2843 >      vy = object->vy - chosenVertex ()->y,
2844 >      vz = object->vz - chosenVertex ()->z,
2845 >      px = object->px,
2846 >      py = object->py,
2847 >      pz = object->pz,
2848 >      pt = object->pt;
2849 >    value = vz - (vx * px + vy * py)/pt * (pz/pt);
2850 >  }
2851 >
2852 >  else if(variable == "deltaRMinSubLeadJet") {
2853 >    // calculate minimum deltaR between track and any other subleading jet  
2854 >    double trkJetDeltaRMin = 99.;  
2855 >    for (uint ijet = 0; ijet<jets->size(); ijet++) {
2856 >      string empty = "";  
2857 >      double isSubLeadingJet = valueLookup(&jets->at(ijet), "disappTrkSubLeadingJetID", "", empty);  
2858 >      if (!isSubLeadingJet) continue;  // only consider jets that pass the subleading jet ID criteria  
2859 >      double jetEta = valueLookup(&jets->at(ijet), "eta", "", empty);
2860 >      double jetPhi = valueLookup(&jets->at(ijet), "phi", "", empty);
2861 >      double trkJetDeltaR = deltaR(object->eta, object->phi, jetEta, jetPhi);  
2862 >      if (trkJetDeltaR < trkJetDeltaRMin) trkJetDeltaRMin = trkJetDeltaR;
2863 >    }
2864 >    value = trkJetDeltaRMin;  
2865 >  }  
2866 >  
2867 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
2868  
2869    else if(variable == "genMatchedPdgId"){
2870      int index = getGenMatchedParticleIndex(object);
# Line 2145 | Line 2885 | OSUAnalysis::valueLookup (const BNtrack*
2885    }
2886    else if(variable == "genMatchedMotherIdReverse"){
2887      int index = getGenMatchedParticleIndex(object);
2888 <    if(index == -1) value = 23;
2889 <    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
2888 >    if(index == -1) value = 24;
2889 >    else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
2890    }
2891    else if(variable == "genMatchedGrandmotherId"){
2892      int index = getGenMatchedParticleIndex(object);
# Line 2158 | Line 2898 | OSUAnalysis::valueLookup (const BNtrack*
2898      }
2899      else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2900    }
2901 +  else if(variable == "genMatchedGrandmotherIdReverse"){
2902 +    int index = getGenMatchedParticleIndex(object);
2903 +    if(index == -1) value = 24;
2904 +    else if(fabs(mcparticles->at(index).motherId) == 15){
2905 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2906 +      if(motherIndex == -1) value = 24;
2907 +      else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2908 +    }
2909 +    else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2910 +  }
2911  
2912  
2913  
2914 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2914 >  else{clog << "WARNING: invalid variable '" << variable << "' for BNtrack collection.\n"; value = -999;}
2915  
2916    value = applyFunction(function, value);
2917  
2918    return value;
2919 < }
2919 > } // end track valueLookup
2920 >
2921  
2922 + //!genjet valueLookup
2923   double
2924 < OSUAnalysis::valueLookup (const BNgenjet* object, string variable, string function){
2924 > OSUAnalysis::valueLookup (const BNgenjet* object, string variable, string function, string &stringValue){
2925  
2926    double value = 0.0;
2927  
# Line 2189 | Line 2941 | OSUAnalysis::valueLookup (const BNgenjet
2941    else if(variable == "charge") value = object->charge;
2942  
2943  
2944 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2944 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2945  
2946    value = applyFunction(function, value);
2947  
2948    return value;
2949   }
2950  
2951 + //!mcparticle valueLookup
2952   double
2953 < OSUAnalysis::valueLookup (const BNmcparticle* object, string variable, string function){
2953 > OSUAnalysis::valueLookup (const BNmcparticle* object, string variable, string function, string &stringValue){
2954  
2955    double value = 0.0;
2956  
# Line 2322 | Line 3075 | OSUAnalysis::valueLookup (const BNmcpart
3075    }
3076  
3077  
3078 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3078 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3079  
3080    value = applyFunction(function, value);
3081  
3082    return value;
3083 < }
3083 > } // end mcparticle valueLookup
3084 >
3085  
3086 + //!primaryvertex valueLookup
3087   double
3088 < OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable, string function){
3088 > OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable, string function, string &stringValue){
3089  
3090    double value = 0.0;
3091  
# Line 2349 | Line 3104 | OSUAnalysis::valueLookup (const BNprimar
3104    else if(variable == "isGood") value = object->isGood;
3105  
3106  
3107 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3107 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3108  
3109    value = applyFunction(function, value);
3110  
3111    return value;
3112   }
3113  
3114 + //!bxlumi valueLookup
3115   double
3116 < OSUAnalysis::valueLookup (const BNbxlumi* object, string variable, string function){
3116 > OSUAnalysis::valueLookup (const BNbxlumi* object, string variable, string function, string &stringValue){
3117  
3118    double value = 0.0;
3119  
# Line 2366 | Line 3122 | OSUAnalysis::valueLookup (const BNbxlumi
3122    else if(variable == "bx_LUMI_now") value = object->bx_LUMI_now;
3123  
3124  
3125 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3125 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3126  
3127    value = applyFunction(function, value);
3128  
3129    return value;
3130   }
3131  
3132 + //!photon valueLookup
3133   double
3134 < OSUAnalysis::valueLookup (const BNphoton* object, string variable, string function){
3134 > OSUAnalysis::valueLookup (const BNphoton* object, string variable, string function, string &stringValue){
3135  
3136    double value = 0.0;
3137  
# Line 2449 | Line 3206 | OSUAnalysis::valueLookup (const BNphoton
3206    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
3207  
3208  
2452
3209  
3210 <  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
3210 >
3211 >  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
3212  
3213    else if(variable == "genMatchedPdgId"){
3214      int index = getGenMatchedParticleIndex(object);
# Line 2473 | Line 3230 | OSUAnalysis::valueLookup (const BNphoton
3230    }
3231    else if(variable == "genMatchedMotherIdReverse"){
3232      int index = getGenMatchedParticleIndex(object);
3233 <    if(index == -1) value = 23;
3234 <    else value = 23 -getPdgIdBinValue(mcparticles->at(index).motherId);
3233 >    if(index == -1) value = 24;
3234 >    else value = 24 -getPdgIdBinValue(mcparticles->at(index).motherId);
3235    }
3236    else if(variable == "genMatchedGrandmotherId"){
3237      int index = getGenMatchedParticleIndex(object);
# Line 2486 | Line 3243 | OSUAnalysis::valueLookup (const BNphoton
3243      }
3244      else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
3245    }
3246 +  else if(variable == "genMatchedGrandmotherIdReverse"){
3247 +    int index = getGenMatchedParticleIndex(object);
3248 +    if(index == -1) value = 24;
3249 +    else if(fabs(mcparticles->at(index).motherId) == 15){
3250 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
3251 +      if(motherIndex == -1) value = 24;
3252 +      else value = 24 - getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
3253 +    }
3254 +    else value = 24 - getPdgIdBinValue(mcparticles->at(index).grandMotherId);
3255 +  }
3256  
3257  
3258 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3258 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3259  
3260    value = applyFunction(function, value);
3261  
3262    return value;
3263 < }
3263 > } // end photon valueLookup
3264 >
3265  
3266 + //!supercluster valueLookup
3267   double
3268 < OSUAnalysis::valueLookup (const BNsupercluster* object, string variable, string function){
3268 > OSUAnalysis::valueLookup (const BNsupercluster* object, string variable, string function, string &stringValue){
3269  
3270    double value = 0.0;
3271  
# Line 2510 | Line 3279 | OSUAnalysis::valueLookup (const BNsuperc
3279    else if(variable == "theta") value = object->theta;
3280  
3281  
3282 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3282 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3283  
3284    value = applyFunction(function, value);
3285  
3286    return value;
3287   }
3288  
3289 + //!trigobj valueLookup
3290 + double
3291 + OSUAnalysis::valueLookup (const BNtrigobj* object, string variable, string function, string &stringValue){
3292 +
3293 +  double value = 0.0;
3294 +
3295 +  if(variable == "pt") value = object->pt;
3296 +  else if(variable == "eta") value = object->eta;
3297 +  else if(variable == "phi") value = object->phi;
3298 +  else if(variable == "px") value = object->px;
3299 +  else if(variable == "py") value = object->py;
3300 +  else if(variable == "pz") value = object->pz;
3301 +  else if(variable == "et") value = object->et;
3302 +  else if(variable == "energy") value = object->energy;
3303 +  else if(variable == "etTotal") value = object->etTotal;
3304 +  else if(variable == "id") value = object->id;
3305 +  else if(variable == "charge") value = object->charge;
3306 +  else if(variable == "isIsolated") value = object->isIsolated;
3307 +  else if(variable == "isMip") value = object->isMip;
3308 +  else if(variable == "isForward") value = object->isForward;
3309 +  else if(variable == "isRPC") value = object->isRPC;
3310 +  else if(variable == "bx") value = object->bx;
3311 +  else if(variable == "filter") {
3312 +    if ((stringValue = object->filter) == "")
3313 +      stringValue = "none";  // stringValue should only be empty if value is filled
3314 +  }
3315 +
3316 +  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3317 +
3318 +  value = applyFunction(function, value);
3319  
3320 +  return value;
3321 + }
3322 +
3323 + //!muon-muon pair valueLookup
3324   double
3325 < OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function){
3325 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function, string &stringValue){
3326  
3327    double value = 0.0;
3328  
# Line 2581 | Line 3384 | OSUAnalysis::valueLookup (const BNmuon*
3384        value = object2->correctedD0;
3385      }
3386  
3387 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3387 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3388 >
3389 >  value = applyFunction(function, value);
3390 >
3391 >  return value;
3392 > } // end muon-muon pair valueLookup
3393 >
3394 >
3395 > //!muon-photon pair valueLookup
3396 > double
3397 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNphoton* object2, string variable, string function, string &stringValue){
3398 >
3399 >  double value = 0.0;
3400 >
3401 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3402 >  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3403 >  else if(variable == "photonEta") value = object2->eta;
3404 >  else if(variable == "photonPt") value = object2->pt;
3405 >  else if(variable == "muonEta") value = object1->eta;
3406 >  else if(variable == "photonPhi") value = object2->phi;
3407 >  else if(variable == "muonPhi") value = object1->phi;
3408 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3409 >  else if(variable == "photonGenMotherId") value = object2->genMotherId;
3410 >  else if(variable == "muonRelPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3411 >  else if(variable == "invMass"){
3412 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3413 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3414 >    value = (fourVector1 + fourVector2).M();
3415 >  }
3416 >  else if(variable == "pt"){
3417 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3418 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3419 >    value = (fourVector1 + fourVector2).Pt();
3420 >  }
3421 >  else if(variable == "threeDAngle")
3422 >    {
3423 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3424 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3425 >      value = (threeVector1.Angle(threeVector2));
3426 >    }
3427 >  else if(variable == "alpha")
3428 >    {
3429 >      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
3430 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3431 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3432 >      value = (pi-threeVector1.Angle(threeVector2));
3433 >    }
3434 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3435 >
3436 >  value = applyFunction(function, value);
3437 >
3438 >  return value;
3439 > }
3440 >
3441 > //!electron-photon pair valueLookup
3442 > double
3443 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNphoton* object2, string variable, string function, string &stringValue){
3444 >
3445 >  double value = 0.0;
3446 >
3447 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3448 >  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3449 >  else if(variable == "photonEta") value = object2->eta;
3450 >  else if(variable == "photonPt") value = object2->pt;
3451 >  else if(variable == "electronEta") value = object1->eta;
3452 >  else if(variable == "photonPhi") value = object2->phi;
3453 >  else if(variable == "electronPhi") value = object1->phi;
3454 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3455 >  else if(variable == "photonGenMotherId") value = object2->genMotherId;
3456 >  else if(variable == "electronRelPFrhoIso") value = ( object1->chargedHadronIsoDR03 + max(0.0, object1->neutralHadronIsoDR03 + object1->photonIsoDR03 - object1->AEffDr03*object1->rhoPrime) ) / object1->pt;
3457 >  else if(variable == "invMass"){
3458 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3459 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3460 >    value = (fourVector1 + fourVector2).M();
3461 >  }
3462 >  else if(variable == "pt"){
3463 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3464 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3465 >    value = (fourVector1 + fourVector2).Pt();
3466 >  }
3467 >  else if(variable == "threeDAngle")
3468 >    {
3469 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3470 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3471 >      value = (threeVector1.Angle(threeVector2));
3472 >    }
3473 >  else if(variable == "alpha")
3474 >    {
3475 >      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
3476 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3477 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3478 >      value = (pi-threeVector1.Angle(threeVector2));
3479 >    }
3480 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3481  
3482    value = applyFunction(function, value);
3483  
3484    return value;
3485   }
3486  
3487 + //!electron-electron pair valueLookup
3488   double
3489 < OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function){
3489 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function, string &stringValue){
3490  
3491    double value = 0.0;
3492  
# Line 2636 | Line 3533 | OSUAnalysis::valueLookup (const BNelectr
3533      value = object2->correctedD0;
3534    }
3535  
3536 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3536 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3537  
3538    value = applyFunction(function, value);
3539  
3540    return value;
3541   }
3542  
3543 + //!electron-muon pair valueLookup
3544   double
3545 < OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function){
3545 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function, string &stringValue){
3546  
3547    double value = 0.0;
3548  
# Line 2702 | Line 3600 | OSUAnalysis::valueLookup (const BNelectr
3600    else if(variable == "muonRelPFdBetaIso"){
3601      value = (object2->pfIsoR04SumChargedHadronPt + max(0.0, object2->pfIsoR04SumNeutralHadronEt + object2->pfIsoR04SumPhotonEt - 0.5*object2->pfIsoR04SumPUPt)) / object2->pt;
3602    }
3603 +  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3604 +  value = applyFunction(function, value);
3605  
3606 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3606 >  return value;
3607 > } // end electron-muon pair valueLookup
3608 >
3609 >
3610 > //!electron-jet pair valueLookup
3611 > double
3612 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNjet* object2, string variable, string function, string &stringValue){
3613  
3614 +  double value = 0.0;
3615 +
3616 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3617 +  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3618 +  else if(variable == "jetEta") value = object2->eta;
3619 +  else if(variable == "jetPhi") value = object2->phi;
3620 +  else if(variable == "electronEta") value = object1->eta;
3621 +  else if(variable == "electronPhi") value = object1->phi;
3622 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3623 +  else if(variable == "invMass"){
3624 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3625 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3626 +    value = (fourVector1 + fourVector2).M();
3627 +  }
3628 +  else if(variable == "pt"){
3629 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3630 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3631 +    value = (fourVector1 + fourVector2).Pt();
3632 +  }
3633 +  else if(variable == "threeDAngle")
3634 +    {
3635 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3636 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3637 +      value = (threeVector1.Angle(threeVector2));
3638 +    }
3639 +  else if(variable == "chargeProduct"){
3640 +    value = object1->charge*object2->charge;
3641 +  }
3642 +
3643 +  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3644    value = applyFunction(function, value);
3645  
3646    return value;
3647   }
3648  
3649 + //!photon-jet pair valueLookup
3650 + double
3651 + OSUAnalysis::valueLookup (const BNphoton* object1, const BNjet* object2, string variable, string function, string &stringValue){
3652 +
3653 +  double value = 0.0;
3654  
3655 < 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));  
3655 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3656    else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3657 <  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);    
3658 <  else if(variable == "invMass"){        
3659 <    fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, electronMass);    
3660 <    fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, electronMass );    
3661 <        
3662 <    value = (fourVector1 + fourVector2).M();    
3657 >  else if(variable == "jetEta") value = object2->eta;
3658 >  else if(variable == "jetPhi") value = object2->phi;
3659 >  else if(variable == "photonEta") value = object1->eta;
3660 >  else if(variable == "photonPhi") value = object1->phi;
3661 >  else if(variable == "photonGenMotherId") value = object1->genMotherId;
3662 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3663 >  else if(variable == "invMass"){
3664 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3665 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3666 >    value = (fourVector1 + fourVector2).M();
3667 >  }
3668 >  else if(variable == "pt"){
3669 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3670 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3671 >    value = (fourVector1 + fourVector2).Pt();
3672    }
3673 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}  
3674 <  value = applyFunction(function, value);        
3675 <  return value;  
3673 >  else if(variable == "threeDAngle")
3674 >    {
3675 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3676 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3677 >      value = (threeVector1.Angle(threeVector2));
3678 >    }
3679 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3680 >  value = applyFunction(function, value);
3681 >
3682 >  return value;
3683   }
3684  
3685 + // track-jet pair valueLookup
3686 + double
3687 + OSUAnalysis::valueLookup (const BNtrack* object1, const BNjet* object2, string variable, string function, string &stringValue){
3688  
3689 +  double value = 0.0;
3690  
3691 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3692 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3693 +
3694 +  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3695 +  value = applyFunction(function, value);
3696 +
3697 +  return value;
3698 +
3699 + }
3700 +
3701 +
3702 +
3703 + // met-jet pair valueLookup
3704   double
3705 < OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrack* object2, string variable, string function){
3705 > OSUAnalysis::valueLookup (const BNmet* object1, const BNjet* object2, string variable, string function, string &stringValue){
3706 >
3707 >  double value = 0.0;
3708 >
3709 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3710 >
3711 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3712 >  value = applyFunction(function, value);
3713 >
3714 >  return value;
3715 >
3716 > }  
3717 >
3718 >
3719 >
3720 > //!muon-jet pair valueLookup
3721 > double
3722 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNjet* object2, string variable, string function, string &stringValue){
3723 >
3724 >  double value = 0.0;
3725 >
3726 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3727 >  else if(variable == "jetEta") value = object2->eta;
3728 >  else if(variable == "relPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3729 >  else if(variable == "jetPt") value = object2->pt;
3730 >  else if(variable == "jetPhi") value = object2->phi;
3731 >  else if(variable == "deltaPt") value = object1->pt - object2->pt;
3732 >  else if(variable == "muonEta") value = object1->eta;
3733 >  else if(variable == "muonPt") value = object1->pt;
3734 >  else if(variable == "muonPhi") value = object1->phi;
3735 >  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);          
3736 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3737 >  else if(variable == "invMass"){
3738 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3739 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3740 >    value = (fourVector1 + fourVector2).M();
3741 >  }
3742 >  else if(variable == "pt"){
3743 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3744 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3745 >    value = (fourVector1 + fourVector2).Pt();
3746 >  }
3747 >  else if(variable == "threeDAngle")
3748 >    {
3749 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3750 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3751 >      value = (threeVector1.Angle(threeVector2));
3752 >    }
3753 >  else if(variable == "chargeProduct"){
3754 >    value = object1->charge*object2->charge;
3755 >  }
3756 >
3757 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3758 >  value = applyFunction(function, value);
3759 >
3760 >  return value;
3761 > }
3762 >
3763 > //!muon-event valueLookup
3764 > double
3765 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNevent* object2, string variable, string function, string &stringValue){
3766 >
3767 >  double value = 0.0;
3768 >
3769 >  if(variable == "muonEta") value = object1->eta;
3770 >  else if(variable == "muonPt") value = object1->pt;
3771 >  else if(variable == "muonPhi") value = object1->phi;
3772 >  else if(variable == "Ht") value = getHt(jets.product());
3773 >  else if(variable == "pthat")   value = object2->pthat;
3774 >  else if(variable == "relPFdBetaIso") value = (object1->pfIsoR04SumChargedHadronPt + max(0.0, object1->pfIsoR04SumNeutralHadronEt + object1->pfIsoR04SumPhotonEt - 0.5*object1->pfIsoR04SumPUPt)) / object1->pt;
3775 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3776 >  value = applyFunction(function, value);
3777 >
3778 >  return value;
3779 > }
3780 > //!jet-jet pair valueLookup
3781 > double
3782 > OSUAnalysis::valueLookup (const BNjet* object1, const BNjet* object2, string variable, string function, string &stringValue){
3783 >
3784 >  double value = 0.0;
3785 >
3786 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3787 >  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3788 >  else if(variable == "deltaPt") value = object1->pt - object2->pt;
3789 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3790 >  else if(variable == "invMass"){
3791 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3792 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3793 >    value = (fourVector1 + fourVector2).M();
3794 >  }
3795 >  else if(variable == "pt"){
3796 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
3797 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
3798 >    value = (fourVector1 + fourVector2).Pt();
3799 >  }
3800 >  else if(variable == "threeDAngle")
3801 >    {
3802 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
3803 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
3804 >      value = (threeVector1.Angle(threeVector2));
3805 >    }
3806 >  else if(variable == "chargeProduct"){
3807 >    value = object1->charge*object2->charge;
3808 >  }
3809 >
3810 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3811 >  value = applyFunction(function, value);
3812 >
3813 >  return value;
3814 > }
3815 >
3816 > //!electron-track pair valueLookup
3817 > double
3818 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3819 >  double electronMass = 0.000511;
3820 >  double value = 0.0;
3821 >  TLorentzVector fourVector1(0, 0, 0, 0);
3822 >  TLorentzVector fourVector2(0, 0, 0, 0);
3823 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3824 >  else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
3825 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3826 >  else if(variable == "invMass"){
3827 >    fourVector1.SetPtEtaPhiM(object1->pt, object1->eta, object1->phi, electronMass);
3828 >    fourVector2.SetPtEtaPhiM(object2->pt, object2->eta, object2->phi, electronMass );
3829 >
3830 >    value = (fourVector1 + fourVector2).M();
3831 >  }
3832 >  else if(variable == "chargeProduct"){
3833 >    value = object1->charge*object2->charge;
3834 >  }
3835 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3836 >  value = applyFunction(function, value);
3837 >  return value;
3838 >
3839 > }
3840 >
3841 >
3842 > //!muon-track pair valueLookup
3843 > double
3844 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3845    double pionMass = 0.140;
3846    double muonMass = 0.106;
3847    double value = 0.0;
# Line 2749 | Line 3856 | OSUAnalysis::valueLookup (const BNmuon*
3856  
3857      value = (fourVector1 + fourVector2).M();
3858    }
3859 +  else if(variable == "chargeProduct"){
3860 +    value = object1->charge*object2->charge;
3861 +  }
3862  
3863 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3863 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3864    value = applyFunction(function, value);
3865    return value;
3866   }
3867  
3868 <
3868 > //!tau-tau pair valueLookup
3869   double
3870 < OSUAnalysis::valueLookup (const BNtau* object1, const BNtau* object2, string variable, string function){
3870 > OSUAnalysis::valueLookup (const BNtau* object1, const BNtau* object2, string variable, string function, string &stringValue){
3871    double value = 0.0;
3872    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3873    else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
# Line 2768 | Line 3878 | OSUAnalysis::valueLookup (const BNtau* o
3878      value = (fourVector1 + fourVector2).M();
3879    }
3880  
3881 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3881 >  else if(variable == "chargeProduct"){
3882 >    value = object1->charge*object2->charge;
3883 >  }
3884 >
3885 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3886    value = applyFunction(function, value);
3887    return value;
3888   }
3889  
3890 <
3890 > //!muon-tau pair valueLookup
3891   double
3892 < OSUAnalysis::valueLookup (const BNmuon* object1, const BNtau* object2, string variable, string function){
3892 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNtau* object2, string variable, string function, string &stringValue){
3893    double value = 0.0;
3894    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3895    else if(variable == "deltaEta") value = fabs(object1->eta - object2->eta);
# Line 2786 | Line 3900 | OSUAnalysis::valueLookup (const BNmuon*
3900      value = (fourVector1 + fourVector2).M();
3901    }
3902  
3903 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3903 >  else if(variable == "chargeProduct"){
3904 >    value = object1->charge*object2->charge;
3905 >  }
3906 >
3907 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3908    value = applyFunction(function, value);
3909    return value;
3910   }
3911  
3912 + //!tau-track pair valueLookup
3913   double
3914 < OSUAnalysis::valueLookup (const BNtau* object1, const BNtrack* object2, string variable, string function){
3914 > OSUAnalysis::valueLookup (const BNtau* object1, const BNtrack* object2, string variable, string function, string &stringValue){
3915    double value = 0.0;
3916    if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
3917    else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3918 +  else if(variable == "chargeProduct"){
3919 +    value = object1->charge*object2->charge;
3920 +  }
3921  
3922 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3922 >  else{clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
3923    value = applyFunction(function, value);
3924    return value;
3925   }
3926  
3927  
3928 <
3928 > //!track-event pair valueLookup
3929   double
3930 < OSUAnalysis::valueLookup (const BNtrack* object1, const BNevent* object2, string variable, string function){
3930 > OSUAnalysis::valueLookup (const BNtrack* object1, const BNevent* object2, string variable, string function, string &stringValue){
3931  
3932    double value = 0.0;
3933    double pMag = sqrt(object1->pt * object1->pt +
3934 <                     object1->pz * object1->pz);  
3934 >                     object1->pz * object1->pz);
3935  
3936    if      (variable == "numPV")                      value = object2->numPV;
3937    else if (variable == "caloTotDeltaRp5")            value =  (object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5);
3938    else if (variable == "caloTotDeltaRp5ByP")         value = ((object1->caloHadDeltaRp5 + object1->caloEMDeltaRp5)/pMag);
3939 <  else if (variable == "caloTotDeltaRp5_RhoCorr")    value = getTrkCaloTotRhoCorr(object1);  
3940 <  else if (variable == "caloTotDeltaRp5ByP_RhoCorr") value = getTrkCaloTotRhoCorr(object1) / pMag;  
3939 >  else if (variable == "caloTotDeltaRp5_RhoCorr")    value = getTrkCaloTotRhoCorr(object1);
3940 >  else if (variable == "caloTotDeltaRp5ByP_RhoCorr") value = getTrkCaloTotRhoCorr(object1) / pMag;
3941  
3942 <  else { std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3942 >  else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3943  
3944    value = applyFunction(function, value);
3945  
3946    return value;
3947  
3948 < }  
3948 > }
3949 >
3950 > //!electron-trigobj pair valueLookup
3951 > double
3952 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNtrigobj* object2, string variable, string function, string &stringValue){
3953 >
3954 >  double value = 0.0;
3955 >
3956 >  if (variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3957 >  else if (variable == "match"){
3958 >    if (deltaR(object1->eta,object1->phi,object2->eta,object2->phi) < 0.2 && abs(object2->id) == 11)
3959 >      stringValue = object2->filter;
3960 >    else
3961 >      stringValue = "none";
3962 >  }
3963 >
3964 >  else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3965 >
3966 >  value = applyFunction(function, value);
3967 >
3968 >  return value;
3969 >
3970 > }
3971 >
3972 > //!muon-trigobj pair valueLookup
3973 > double
3974 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNtrigobj* object2, string variable, string function, string &stringValue){
3975 >
3976 >  double value = 0.0;
3977 >
3978 >  if (variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
3979 >
3980 >  else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
3981 >
3982 >  value = applyFunction(function, value);
3983 >
3984 >  return value;
3985 >
3986 > }
3987 >
3988 > //!stop valueLookup
3989 > double
3990 > OSUAnalysis::valueLookup (const BNstop* object, string variable, string function, string &stringValue){
3991 >
3992 >
3993 >  double value = 0.0;
3994 >
3995 >  if(variable == "ctau") value = object->ctau;
3996 >
3997 >  else if (variable == "d0"){
3998 >    double vx = object->vx - chosenVertex ()->x,
3999 >      vy = object->vy - chosenVertex ()->y,
4000 >      px = object->px,
4001 >      py = object->py,
4002 >      pt = object->pt;
4003 >    value = (-vx * py + vy * px) / pt;
4004 >  }
4005 >
4006 >  else if (variable == "dz"){
4007 >    double vx = object->vx - chosenVertex ()->x,
4008 >      vy = object->vy - chosenVertex ()->y,
4009 >      vz = object->vz - chosenVertex ()->z,
4010 >      px = object->px,
4011 >      py = object->py,
4012 >      pz = object->pz,
4013 >      pt = object->pt;
4014 >    value = vz - (vx * px + vy * py)/pt * (pz/pt);
4015 >  }
4016 >
4017 >  else if (variable == "minD0"){
4018 >    double minD0=999;
4019 >    for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4020 >      double vx = object->vx - vertex->x,
4021 >        vy = object->vy - vertex->y,
4022 >        px = object->px,
4023 >        py = object->py,
4024 >        pt = object->pt;
4025 >      value = (-vx * py + vy * px) / pt;
4026 >      if(abs(value) < abs(minD0)) minD0 = value;
4027 >    }
4028 >    value = minD0;
4029 >  }
4030 >  else if (variable == "minDz"){
4031 >    double minDz=999;
4032 >    for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4033 >      double vx = object->vx - vertex->x,
4034 >        vy = object->vy - vertex->y,
4035 >        vz = object->vz - vertex->z,
4036 >        px = object->px,
4037 >        py = object->py,
4038 >        pz = object->pz,
4039 >        pt = object->pt;
4040 >      value = vz - (vx * px + vy * py)/pt * (pz/pt);
4041 >      if(abs(value) < abs(minDz)) minDz = value;
4042 >    }
4043 >    value = minDz;
4044 >  }
4045 >  else if(variable == "distToVertex"){
4046 >    value = sqrt((object->vx-chosenVertex()->x)*(object->vx-chosenVertex()->x) + \
4047 >                 (object->vy-chosenVertex()->y)*(object->vy-chosenVertex()->y) + \
4048 >                 (object->vz-chosenVertex()->z)*(object->vz-chosenVertex()->z));
4049 >  }
4050 >  else if (variable == "minDistToVertex"){
4051 >    double minDistToVertex=999;
4052 >    for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4053 >      value = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4054 >                   (object->vy-vertex->y)*(object->vy-vertex->y) + \
4055 >                   (object->vz-vertex->z)*(object->vz-vertex->z));
4056 >
4057 >      if(abs(value) < abs(minDistToVertex)) minDistToVertex = value;
4058 >    }
4059 >    value = minDistToVertex;
4060 >  }
4061 >  else if (variable == "distToVertexDifference"){
4062 >    double minDistToVertex=999;
4063 >    for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4064 >      value = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4065 >                   (object->vy-vertex->y)*(object->vy-vertex->y) + \
4066 >                   (object->vz-vertex->z)*(object->vz-vertex->z));
4067 >
4068 >      if(abs(value) < abs(minDistToVertex)) minDistToVertex = value;
4069 >    }
4070 >    double distToChosenVertex = sqrt((object->vx-chosenVertex()->x)*(object->vx-chosenVertex()->x) + \
4071 >                                     (object->vy-chosenVertex()->y)*(object->vy-chosenVertex()->y) + \
4072 >                                     (object->vz-chosenVertex()->z)*(object->vz-chosenVertex()->z));
4073 >
4074 >    value = distToChosenVertex - minDistToVertex;
4075 >  }
4076 >
4077 >  else if (variable == "closestVertexRank"){
4078 >    double minDistToVertex=999;
4079 >    int vertex_rank = 0;
4080 >    for(BNprimaryvertexCollection::const_iterator vertex = primaryvertexs->begin (); vertex != primaryvertexs->end (); vertex++){
4081 >      vertex_rank++;
4082 >      int dist = sqrt((object->vx-vertex->x)*(object->vx-vertex->x) + \
4083 >                      (object->vy-vertex->y)*(object->vy-vertex->y) + \
4084 >                      (object->vz-vertex->z)*(object->vz-vertex->z));
4085 >
4086 >      if(abs(dist) < abs(minDistToVertex)){
4087 >        value = vertex_rank;
4088 >        minDistToVertex = dist;
4089 >      }
4090 >    }
4091 >  }
4092 >
4093 >  else if (variable == "decaysToTau"){
4094 >    value = abs (object->daughter0Id) == 15 || abs (object->daughter1Id) == 15;
4095 >  }
4096 >
4097 >
4098 >
4099 >
4100 >  else { clog << "WARNING: invalid variable '" << variable << "'\n"; value = -999; }
4101 >
4102 >  value = applyFunction(function, value);
4103 >
4104 >  return value;
4105 >
4106 > } // end stop valueLookup
4107 >
4108 >
4109 >
4110 >
4111  
4112   // Calculate the number of tracks in cone of DeltaR<0.5 around track1.
4113   // Return true iff no other tracks are found in this cone.
# Line 2838 | Line 4122 | OSUAnalysis::getTrkIsIso (const BNtrack*
4122  
4123   }
4124  
4125 + //calculate the scalar sum of Jet Pt in the event.
4126 + double
4127 + OSUAnalysis::getHt (const BNjetCollection* jetColl){
4128 +  double Ht = 0;
4129 +  for(BNjetCollection::const_iterator jet = jetColl->begin(); jet !=jetColl->end(); jet++){
4130 +    Ht += abs(jet->pt);
4131 +  }
4132 +  return Ht;
4133 + }
4134  
4135   double
4136   OSUAnalysis::getTrkPtRes (const BNtrack* track1){
# Line 2872 | Line 4165 | OSUAnalysis::getTrkPtTrue (const BNtrack
4165  
4166   double
4167   OSUAnalysis::getTrkCaloTotRhoCorr(const BNtrack* track) {
4168 <  // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.  
4169 <  if (!useTrackCaloRhoCorr_) return -99;  
4170 <  // if (!rhokt6CaloJetsHandle_) {
4171 <  //   cout << "ERROR [getTrkCaloTotRhoCorr]:  The collection rhokt6CaloJetsHandle is not available!" << endl;  
4172 <  //   return -99;  
4168 >  // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.
4169 >  if (!useTrackCaloRhoCorr_) return -99;
4170 >  // if (!rhokt6CaloJetsHandle_) {
4171 >  //   clog << "ERROR [getTrkCaloTotRhoCorr]:  The collection rhokt6CaloJetsHandle is not available!" << endl;
4172 >  //   return -99;
4173    // }
4174 <  double radDeltaRCone = 0.5;  
4175 <  double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);  // Define effective area as pi*r^2, where r is radius of DeltaR cone.  
4176 <  double rawCaloTot = track->caloHadDeltaRp5 + track->caloEMDeltaRp5;  
4177 <  double caloTotRhoCorrCalo = TMath::Max(0., rawCaloTot - rhoCorr_kt6CaloJets);  
4178 <  return caloTotRhoCorrCalo;  
4179 <
4174 >  double radDeltaRCone = 0.5;
4175 >  double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);  // Define effective area as pi*r^2, where r is radius of DeltaR cone.
4176 >  double rawCaloTot = track->caloHadDeltaRp5 + track->caloEMDeltaRp5;
4177 >  double caloTotRhoCorrCalo = TMath::Max(0., rawCaloTot - rhoCorr_kt6CaloJets);
4178 >  return caloTotRhoCorrCalo;
4179 >
4180 > }
4181 >
4182 > double
4183 > OSUAnalysis::getTrkDepTrkRp5RhoCorr(const BNtrack* track) {
4184 >  // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.              
4185 >  if (!useTrackCaloRhoCorr_) return -99;
4186 >  // if (!rhokt6CaloJetsHandle_) {                                                                                                    
4187 >  //   clog << "ERROR [getTrkCaloTotRhoCorr]:  The collection rhokt6CaloJetsHandle is not available!" << endl;                        
4188 >  //   return -99;                                                                                                                    
4189 >  // }                                                                                                                                
4190 >  double radDeltaRCone = 0.5;
4191 >  double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);  // Define effective area as pi*r^2, where r is radius of DeltaR cone.                                                                                                          
4192 >  double rawDepTrkRp5 = track->depTrkRp5;
4193 > double depTrkRp5RhoCorr = TMath::Max(0., rawDepTrkRp5 - rhoCorr_kt6CaloJets);
4194 > return depTrkRp5RhoCorr;
4195 >
4196 > }
4197 >
4198 > double
4199 > OSUAnalysis::getTrkDepTrkRp3RhoCorr(const BNtrack* track) {
4200 >  // Return the pile-up (rho) corrected isolation energy, i.e., the total calorimeter energy around the candidate track.                                                
4201 >  if (!useTrackCaloRhoCorr_) return -99;
4202 >  // if (!rhokt6CaloJetsHandle_) {                                                                                                                                      
4203 >  //   clog << "ERROR [getTrkCaloTotRhoCorr]:  The collection rhokt6CaloJetsHandle is not available!" << endl;                                                          
4204 >  //   return -99;                                                                                                                                                      
4205 >  // }                                                                                                                                                                  
4206 >  double radDeltaRCone = 0.3;
4207 > // Define effective area as pi*r^2, where r is radius of DeltaR cone
4208 >  double rhoCorr_kt6CaloJets = *rhokt6CaloJetsHandle_ * TMath::Pi() * pow(radDeltaRCone, 2);  
4209 >  double rawDepTrkRp3 = track->depTrkRp3;
4210 >  double depTrkRp3RhoCorr = TMath::Max(0., rawDepTrkRp3 - rhoCorr_kt6CaloJets);
4211 >  return depTrkRp3RhoCorr;
4212 >  
4213   }
4214  
4215  
# Line 2897 | Line 4223 | OSUAnalysis::WriteDeadEcal (){
4223    double etaEcal, phiEcal;
4224    ifstream DeadEcalFile(deadEcalFile_);
4225    if(!DeadEcalFile) {
4226 <    cout << "Error: DeadEcalFile has not been found." << endl;
4226 >    clog << "Error: DeadEcalFile has not been found." << endl;
4227      return;
4228    }
4229    if(DeadEcalVec.size()!= 0){
4230 <    cout << "Error: DeadEcalVec has a nonzero size" << endl;
4230 >    clog << "Error: DeadEcalVec has a nonzero size" << endl;
4231      return;
4232    }
4233    while(!DeadEcalFile.eof())
# Line 2912 | Line 4238 | OSUAnalysis::WriteDeadEcal (){
4238        newChan.phiEcal = phiEcal;
4239        DeadEcalVec.push_back(newChan);
4240      }
4241 <  if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl;
4241 >  if(DeadEcalVec.size() == 0) clog << "Warning: No dead Ecal channels have been found." << endl;
4242   }
4243  
4244   //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
# Line 2921 | Line 4247 | OSUAnalysis::getTrkIsMatchedDeadEcal (co
4247    double deltaRLowest = 999;
4248    int value = 0;
4249    if (DeadEcalVec.size() == 0) WriteDeadEcal();
4250 <  for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
4250 >  for(vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
4251      double eta = ecal->etaEcal;
4252      double phi = ecal->phiEcal;
4253 <    double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
4253 >    double deltaRtemp = deltaR(eta, phi, track1->eta, track1->phi);
4254      if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
4255    }
4256    if (deltaRLowest<0.05) {value = 1;}
# Line 2932 | Line 4258 | OSUAnalysis::getTrkIsMatchedDeadEcal (co
4258    return value;
4259   }
4260  
4261 < // Returns the smallest DeltaR between the object and any generated true particle in the event.  
4261 > // Returns the smallest DeltaR between the object and any generated true particle in the event.
4262   template <class InputObject>
4263   double OSUAnalysis::getGenDeltaRLowest(InputObject object){
4264    double genDeltaRLowest = 999.;
# Line 2952 | Line 4278 | OSUAnalysis::applyFunction(string functi
4278    else if(function == "log") value = log10(value);
4279  
4280    else if(function == "") value = value;
4281 <  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
4281 >  else{clog << "WARNING: invalid function '" << function << "'\n";}
4282  
4283    return value;
4284  
# Line 2961 | Line 4287 | OSUAnalysis::applyFunction(string functi
4287  
4288   template <class InputCollection>
4289   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
4290 <
4291 <  if (currentCut.inputCollection.find("pair")!=std::string::npos)  {
4292 <    string obj1, obj2;
4290 >
4291 >  if (verbose_>2) clog << "  Beginning setObjectFlags for cut " << currentCutIndex << ": " << currentCut.name
4292 >                       << ", inputType=" << inputType
4293 >                       << endl;  
4294 >  if (currentCut.inputCollection.find("pair")!=string::npos)  {
4295 >    string obj1, obj2;
4296      getTwoObjs(currentCut.inputCollection, obj1, obj2);
4297      if (inputType==obj1 ||
4298 <        inputType==obj2) {
4298 >        inputType==obj2) {
4299        // Do not add a cut to individualFlags or cumulativeFlags, if the cut is on a paired collection,
4300 <      // and the inputType is a member of the pair.  
4301 <      // The cut will instead be applied when the setObjectFlags() is called for the paired collection.  
4302 <      // For example, if currentCut.inputCollection==electron-muon pairs,
4303 <      // then the flags should not be set here when inputType==muons or inputType==electrons.  
4304 <      return;
4305 <    }  
4306 <  }    
4300 >      // and the inputType is a member of the pair.
4301 >      // The cut will instead be applied when the setObjectFlags() is called for the paired collection.
4302 >      // For example, if currentCut.inputCollection==electron-muon pairs,
4303 >      // then the flags should not be set here when inputType==muons or inputType==electrons.
4304 >      return;
4305 >    }
4306 >  }
4307  
4308    for (uint object = 0; object != inputCollection->size(); object++){
4309  
4310 <    bool decision = true;//object passes if this cut doesn't cut on that type of object
4310 >    bool cutDecision = true;//object passes if this cut doesn't cut on that type of object
4311 >    bool plotDecision = true;
4312  
4313      if(currentCut.inputCollection == inputType){
4314  
4315        vector<bool> subcutDecisions;
4316        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
4317 <        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
4318 <        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4317 >        string stringValue = "";
4318 >        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex), stringValue);
4319 >        if (stringValue == "") subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4320 >        else subcutDecisions.push_back(evaluateComparison(stringValue,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutStringValues.at(subcutIndex)));
4321  
4322        }
4323 <      if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
4323 >      if(currentCut.numSubcuts == 1) cutDecision = subcutDecisions.at(0);
4324        else{
4325          bool tempDecision = true;
4326          for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
# Line 2997 | Line 4329 | void OSUAnalysis::setObjectFlags(cut &cu
4329            else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
4330              tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
4331          }
4332 <        decision = tempDecision;
4332 >        cutDecision = tempDecision;
4333        }
4334 +      //invert the cut if this cut is a veto
4335 +      if(currentCut.isVeto) cutDecision = !cutDecision;
4336 +      plotDecision = cutDecision;
4337      }
4338  
4339 <    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
4339 >    individualFlags.at(inputType).at(currentCutIndex).push_back(make_pair(cutDecision,plotDecision));
4340  
4341 <    //set flags for objects that pass each cut AND all the previous cuts
4342 <    bool previousCumulativeFlag = true;
4341 >    //set flags for objects that pass this cut AND all the previous cuts
4342 >    bool previousCumulativeCutFlag = true;
4343      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4344 <      if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(object)) previousCumulativeFlag = true;
4345 <      else{ previousCumulativeFlag = false; break;}
4344 >      if(previousCumulativeCutFlag && individualFlags.at(inputType).at(previousCutIndex).at(object).first) previousCumulativeCutFlag = true;
4345 >      else{ previousCumulativeCutFlag = false; break;}
4346      }
4347 <    cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
4347 >    previousCumulativeCutFlag = previousCumulativeCutFlag && cutDecision;
4348 >    bool previousCumulativePlotFlag = true;
4349 >    for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4350 >      if(previousCumulativePlotFlag && individualFlags.at(inputType).at(previousCutIndex).at(object).second) previousCumulativePlotFlag = true;
4351 >      else{ previousCumulativePlotFlag = false; break;}
4352 >    }
4353 >    previousCumulativePlotFlag = previousCumulativePlotFlag && plotDecision;
4354 >
4355 >    cumulativeFlags.at(inputType).at(currentCutIndex).push_back(make_pair(previousCumulativeCutFlag,previousCumulativePlotFlag));
4356  
4357    }
4358  
4359 < }
4359 > } // end void OSUAnalysis::setObjectFlags
4360  
4361  
4362   template <class InputCollection1, class InputCollection2>
4363 < void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
4364 <                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
4365 <
4363 > void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags,
4364 >                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, string inputType){
4365 >  // This function sets the flags for the paired object collection.  
4366 >  // If the cut is applying on the given paired object collection, then the flags for the single object collections are also set.
4367 >  // If not, then the flags for the paired object collection are taken as the AND of the flags for each single object collection.  
4368 >
4369 >  if (verbose_>2) clog << "  Beginning setObjectFlags for cut=" << currentCut.name
4370 >                       << ", inputType=" << inputType
4371 >                       << endl;  
4372  
4373    bool sameObjects = false;
4374 <  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
4374 >  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;  // FIXME:  is sameObjects just the not of isTwoTypesOfObject?  If so, it's redundant.  
4375  
4376 <  // Get the strings for the two objects that make up the pair.  
4377 <  string obj1Type, obj2Type;
4376 >  // Get the strings for the two objects that make up the pair.
4377 >  string obj1Type, obj2Type;
4378    getTwoObjs(inputType, obj1Type, obj2Type);
4379    bool isTwoTypesOfObject = true;
4380 <  if (obj1Type==obj2Type) isTwoTypesOfObject = false;  
4380 >  if (obj1Type==obj2Type) isTwoTypesOfObject = false;
4381  
4382 <  // Initialize the flags for individual objects to all be false, if the cut is on the pair.  
4383 <  // Set them to true later, if any paired object passes (in which case both of its constituents should pass).  
4384 <  if (currentCut.inputCollection == inputType) {    
4382 >  // Initialize the flags for individual objects to all be false, if the cut is on the pair.
4383 >  // Set them to true later, if any paired object passes (in which case both of its constituents should pass).
4384 >  if (currentCut.inputCollection == inputType) {
4385      for (uint object1 = 0; object1 != inputCollection1->size(); object1++) {
4386 <      individualFlags.at(obj1Type).at(currentCutIndex).push_back(false);
4387 <      cumulativeFlags.at(obj1Type).at(currentCutIndex).push_back(false);
4386 >      individualFlags.at(obj1Type).at(currentCutIndex).push_back(make_pair(false,false));
4387 >      cumulativeFlags.at(obj1Type).at(currentCutIndex).push_back(make_pair(false,false));
4388      }
4389 <    if (isTwoTypesOfObject) { // Only initialize the second object if it is different from the first.  
4389 >    if (isTwoTypesOfObject) { // Only initialize the second object if it is different from the first.
4390        for (uint object2 = 0; object2 != inputCollection2->size(); object2++)  {
4391 <        individualFlags.at(obj2Type).at(currentCutIndex).push_back(false);
4392 <        cumulativeFlags.at(obj2Type).at(currentCutIndex).push_back(false);
4391 >        individualFlags.at(obj2Type).at(currentCutIndex).push_back(make_pair(false,false));
4392 >        cumulativeFlags.at(obj2Type).at(currentCutIndex).push_back(make_pair(false,false));
4393        }
4394      }
4395    }
4396 <  
4396 >
4397    int counter = 0;
4398  
4399    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
# Line 3053 | Line 4402 | void OSUAnalysis::setObjectFlags(cut &cu
4402        if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
4403  
4404  
4405 <      bool decision = true;//object passes if this cut doesn't cut on that type of object
4405 >      bool cutDecision  = true;//object passes if this cut doesn't cut on that type of object
4406 >      bool plotDecision = true;
4407  
4408 +      // Determine whether each pair passes the cut, only if inputCollection is the same as the inputType.  
4409        if(currentCut.inputCollection == inputType){
4410  
4411          vector<bool> subcutDecisions;
4412          for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
4413 <          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
4414 <          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4413 >          string stringValue = "";
4414 >          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex), stringValue);
4415 >          if (verbose_>1) clog << currentCut.variables.at(subcutIndex) << " = " << value
4416 >                               << endl;  
4417 >          if (stringValue == "") subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
4418 >          else subcutDecisions.push_back(evaluateComparison(stringValue,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutStringValues.at(subcutIndex)));
4419          }
4420  
4421 <        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
4421 >        if(currentCut.numSubcuts == 1) cutDecision = subcutDecisions.at(0);
4422          else{
4423            bool tempDecision = subcutDecisions.at(0);
4424            for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
# Line 3072 | Line 4427 | void OSUAnalysis::setObjectFlags(cut &cu
4427              else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
4428                tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
4429            }
4430 <          decision = tempDecision;
4430 >          cutDecision = tempDecision;
4431          }
4432 +        //invert the cut if this cut is a veto
4433 +        if (currentCut.isVeto) cutDecision = !cutDecision;
4434 +        plotDecision = cutDecision;
4435 +
4436 +        if (verbose_>1) clog << " cutDecision = " << cutDecision
4437 +                             << "; for currentCut.inputCollection = " << currentCut.inputCollection
4438 +                             << "; object1 (" << obj1Type << ") = " << object1
4439 +                             << "; object2 (" << obj2Type << ") = " << object2
4440 +                             << endl;  
4441 +        
4442 +        if (cutDecision) {  // only set the flags for the individual objects if the pair object is being cut on
4443 +          individualFlags.at(obj1Type).at(currentCutIndex).at(object1).first = true;
4444 +          individualFlags.at(obj2Type).at(currentCutIndex).at(object2).first = true;
4445 +        }
4446 +        if (plotDecision) {  // only set the flags for the individual objects if the pair object is being cut on
4447 +          individualFlags.at(obj1Type).at(currentCutIndex).at(object1).second = true;
4448 +          individualFlags.at(obj2Type).at(currentCutIndex).at(object2).second = true;
4449 +        }
4450 +        
4451 +
4452 +      } // if(currentCut.inputCollection == inputType){
4453 +
4454 +      // The individualFlags will be true if the inputCollection is not the same as the inputType.
4455 +      // They are also independent of the previous flags on the single objects.  
4456 +      individualFlags.at(inputType).at(currentCutIndex).push_back(make_pair(cutDecision,plotDecision));  
4457 +
4458 +
4459 +
4460 +      // ************************************
4461 +      // Determine cumulative flags
4462 +      // ************************************
4463 +      // determine whether this paired object passes this cut AND all previous cuts
4464 +      bool previousCumulativeCutFlag = true;
4465 +      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4466 +        if(previousCumulativeCutFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter).first) previousCumulativeCutFlag = true;
4467 +        else{ previousCumulativeCutFlag = false; break;}
4468        }
4469 <      // 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 <      }  
4469 >      previousCumulativeCutFlag = previousCumulativeCutFlag && cutDecision;
4470  
4471 <      //set flags for objects that pass each cut AND all the previous cuts
3087 <      bool previousCumulativeFlag = true;
4471 >      bool previousCumulativePlotFlag = true;
4472        for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
4473 <        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
4474 <        else{ previousCumulativeFlag = false; break;}
4473 >        if(previousCumulativePlotFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter).second) previousCumulativePlotFlag = true;
4474 >        else{ previousCumulativePlotFlag = false; break;}
4475        }
4476 <      //apply flags for the components of the composite object as well
4477 <      bool currentCumulativeFlag = true;
4478 <      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
4479 <      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
4480 <      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
4481 <      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
4482 <      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
4483 <      if (currentCumulativeFlag && currentCut.inputCollection == inputType) {  // only set the flags for the individual objects if the pair object is being cut on  
4484 <        cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1) = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1);  
4485 <        cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2) = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2);  
4476 >      previousCumulativePlotFlag = previousCumulativePlotFlag && plotDecision;
4477 >      
4478 >      // Get the index for the flags of each of the single objects in the pair.  Usually this is the index of the previous cut, i.e., currentCutIndex-1.    
4479 >      int cutIdxFlagsObj1 = max(int(currentCutIndex-1),0);
4480 >      int cutIdxFlagsObj2 = max(int(currentCutIndex-1),0);
4481 >      // If the inputCollection of the cut is not equal to the inputType but the inputCollection includes objects that are contained in inputType, then use the currentCutIndex for those collections.
4482 >      // For example, if the inputType is jet-jet pairs, and the inputCollection is track-jet pairs, then use currentCutIndex for cutIdxFlagsObj{1,2}, i.e., for both jets.  
4483 >      // For example, if the inputType is jets, and the inputCollection is track-jet pairs, then use currentCutIndex for cutIdxFlagsObj2 (jets) and currentCutIndex-1 for cutIdxFlagsObj1 (tracks).  
4484 >      if (currentCut.inputCollection != inputType) {  
4485 >        if (currentCut.inputCollection.find(obj1Type)!=string::npos) cutIdxFlagsObj1 = currentCutIndex;
4486 >        if (currentCut.inputCollection.find(obj2Type)!=string::npos) cutIdxFlagsObj2 = currentCutIndex;
4487 >      }  
4488 >      flagPair flags1 = cumulativeFlags.at(obj1Type).at(cutIdxFlagsObj1);  // flag for input collection 1
4489 >      flagPair flags2 = cumulativeFlags.at(obj2Type).at(cutIdxFlagsObj2);  // flag for input collection 2
4490 >      
4491 >      // The cumulative flag is only true if the paired object cumulative flag is true, and if the single object cumulative flags are true.  
4492 >      bool currentCumulativeCutFlag = true;
4493 >      bool currentCumulativePlotFlag = true;
4494 >      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag;
4495 >      else if(flags1.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag && flags2.at(object2).first;
4496 >      else if(flags2.size() == 0) currentCumulativeCutFlag = previousCumulativeCutFlag && flags1.at(object1).first;
4497 >      else currentCumulativeCutFlag = previousCumulativeCutFlag && flags1.at(object1).first && flags2.at(object2).first;
4498 >
4499 >      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag;
4500 >      else if(flags1.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag && flags2.at(object2).second;
4501 >      else if(flags2.size() == 0) currentCumulativePlotFlag = previousCumulativePlotFlag && flags1.at(object1).second;
4502 >      else currentCumulativePlotFlag = previousCumulativePlotFlag && flags1.at(object1).first && flags2.at(object2).second;
4503 >
4504 >      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(make_pair(currentCumulativeCutFlag,currentCumulativePlotFlag));  // Set the flag for the paired object
4505 >
4506 >
4507 >      if (currentCumulativeCutFlag && currentCut.inputCollection == inputType) {  // Set the flags for the individual objects if the paired object is being cut on.  
4508 >        cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1).first  = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1, "cut");
4509 >        cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2).first  = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "cut");
4510 >        cumulativeFlags.at(obj1Type).at(currentCutIndex).at(object1).second = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj1Type, object1, "plot");
4511 >        cumulativeFlags.at(obj2Type).at(currentCutIndex).at(object2).second = true && getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "plot");
4512 >
4513 >        if (verbose_>1) clog << " previousCumulativeCutFlag for object1 = " << getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object1, "cut") << endl;  
4514 >        if (verbose_>1) clog << " previousCumulativeCutFlag for object2 = " << getPreviousCumulativeFlags(currentCutIndex, individualFlags, obj2Type, object2, "cut") << endl;  
4515 >
4516        }
4517 +
4518        counter++;
4519  
4520 <    } // end   for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
4521 <  }  // end   for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
4520 >    } // end   for (uint object2 = 0; object2 != inputCollection2->size(); object2++)
4521 >  }  // end   for (uint object1 = 0; object1 != inputCollection1->size(); object1++)
4522 >
4523 > } // end void OSUAnalysis::setObjectFlags
4524  
3108 }
4525  
4526 + bool OSUAnalysis::getPreviousCumulativeFlags(uint currentCutIndex, flagMap &individualFlags, string obj1Type, uint object1, string flagType) {
4527 +  // Return true iff for the collection obj1Type, the element with index object1 has individal flags set to true for
4528 +  // all cuts up to currentCutIndex
4529 +  bool previousCumulativeFlag = true;
4530 +  for (uint previousCutIndex = 0; previousCutIndex < currentCutIndex; previousCutIndex++) {
4531 +    bool tempFlag = false;
4532 +    if(flagType == "cut") tempFlag = individualFlags.at(obj1Type).at(previousCutIndex).at(object1).first;
4533 +    else if(flagType == "plot") tempFlag = individualFlags.at(obj1Type).at(previousCutIndex).at(object1).second;
4534  
4535 < bool OSUAnalysis::getPreviousCumulativeFlags(uint currentCutIndex, flagMap &individualFlags, string obj1Type, uint object1) {
4536 <  // Return true iff for the collection obj1Type, the element with index object1 has individal flags set to true for
4537 <  // 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;  
4535 >    if (previousCumulativeFlag && tempFlag) previousCumulativeFlag = true;
4536 >    else {
4537 >      previousCumulativeFlag = false; break;
4538      }
4539    }
4540 <  return previousCumulativeFlag;  
4541 < }  
4540 >  return previousCumulativeFlag;
4541 > }
4542 >
4543 >
4544 >
4545 >
4546 > template <class InputCollection>
4547 > void OSUAnalysis::assignTreeBranch(BranchSpecs parameters, InputCollection inputCollection, flagPair flags){
4548 >  // This function is similar to fill1DHistogram(), but instead of filling a histogram it assigns a value to a variable for the BNTree
4549 >
4550 >  if (BNTreeBranchVals_.count(parameters.name)==0) clog << "Error[assignTreeBranch]:  trying to assign value to " << parameters.name << " that does not have a branch set up.  Will likely seg fault." << endl;
4551 >  for (uint object = 0; object != inputCollection->size(); object++) {
4552 >
4553 >    if (!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4554 >
4555 >    string inputVariable = parameters.inputVariable;
4556 >    string function = "";
4557 >    string stringValue = "";
4558 >    double value = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4559 >    BNTreeBranchVals_.at(parameters.name).push_back(value);
4560 >
4561 >  }
4562 > }
4563  
4564  
4565   template <class InputCollection>
4566 < void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
4566 > void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection, flagPair flags, double scaleFactor){
4567 >
4568 >  if (verbose_>2) clog << "  Filling histogram for " << parameters.name << endl;  
4569  
4570    for (uint object = 0; object != inputCollection->size(); object++){
4571  
4572 <    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
4572 >    if(!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4573  
4574      string currentString = parameters.inputVariables.at(0);
4575      string inputVariable = "";
4576      string function = "";
4577 <    if(currentString.find("(")==std::string::npos){
4577 >    if(currentString.find("(")==string::npos){
4578        inputVariable = currentString;// variable to cut on
4579      }
4580      else{
# Line 3141 | Line 4583 | void OSUAnalysis::fill1DHistogram(TH1* h
4583        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4584      }
4585  
4586 <    double value = valueLookup(&inputCollection->at(object), inputVariable, function);
4586 >    string stringValue = "";
4587 >    double value = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4588      histo->Fill(value,scaleFactor);
4589  
4590      if (printEventInfo_) {
4591 <      // Write information about event to screen, for testing purposes.  
4592 <      cout << "  Info for event:  value for histogram " << histo->GetName() << ":  " << value << endl;  
4591 >      // Write information about event to screen, for testing purposes.
4592 >      clog << "  Info for event:  value for histogram " << histo->GetName() << ":  " << value << " (object number " << object << ")" << endl;
4593      }
4594 <    
4594 >
4595    }
4596   }
4597  
4598   template <class InputCollection1, class InputCollection2>
4599 < void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
4599 > void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2,
4600 >                                  flagPair pairFlags, double scaleFactor){
4601  
4602    bool sameObjects = false;
4603    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
# Line 3165 | Line 4609 | void OSUAnalysis::fill1DHistogram(TH1* h
4609        if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
4610  
4611        pairCounter++;
4612 <      //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;
3171 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
4612 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter).second) continue;
4613  
4614        string currentString = parameters.inputVariables.at(0);
4615        string inputVariable = "";
4616        string function = "";
4617 <      if(currentString.find("(")==std::string::npos){
4617 >      if(currentString.find("(")==string::npos){
4618          inputVariable = currentString;// variable to cut on
4619        }
4620        else{
# Line 3182 | Line 4623 | void OSUAnalysis::fill1DHistogram(TH1* h
4623          inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4624        }
4625  
4626 <      double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
4626 >      string stringValue = "";
4627 >      double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4628        histo->Fill(value,scaleFactor);
4629  
4630 +      if (printEventInfo_) {
4631 +        // Write information about event to screen, for testing purposes.
4632 +        clog << "  Info for event:  value for histogram " << histo->GetName() << ":  " << value
4633 +             << " (object1 number " << object1 << "), "
4634 +             << " (object2 number " << object2 << ")"
4635 +             << endl;
4636 +      }
4637 +      
4638      }
4639    }
4640  
# Line 3192 | Line 4642 | void OSUAnalysis::fill1DHistogram(TH1* h
4642  
4643  
4644   template <class InputCollection>
4645 < void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
4645 > void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection, flagPair flags, double scaleFactor){
4646  
4647    for (uint object = 0; object != inputCollection->size(); object++){
4648  
4649 <    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
4649 >    if(!plotAllObjectsInPassingEvents_ && !flags.at(object).second) continue;
4650  
4651 +    string stringValue = "";
4652      string currentString = parameters.inputVariables.at(0);
4653      string inputVariable = "";
4654      string function = "";
4655 <    if(currentString.find("(")==std::string::npos){
4655 >    if(currentString.find("(")==string::npos){
4656        inputVariable = currentString;// variable to cut on
4657      }
4658      else{
# Line 3209 | Line 4660 | void OSUAnalysis::fill2DHistogram(TH2* h
4660        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4661        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4662      }
4663 <    double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
4663 >    double valueX = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4664  
4665      currentString = parameters.inputVariables.at(1);
4666      inputVariable = "";
4667      function = "";
4668 <    if(currentString.find("(")==std::string::npos){
4668 >    if(currentString.find("(")==string::npos){
4669        inputVariable = currentString;// variable to cut on
4670      }
4671      else{
# Line 3223 | Line 4674 | void OSUAnalysis::fill2DHistogram(TH2* h
4674        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4675      }
4676  
4677 <    double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
4677 >    double valueY = valueLookup(&inputCollection->at(object), inputVariable, function, stringValue);
4678  
4679      histo->Fill(valueX,valueY,scaleFactor);
4680  
# Line 3232 | Line 4683 | void OSUAnalysis::fill2DHistogram(TH2* h
4683   }
4684  
4685   template <class InputCollection1, class InputCollection2>
4686 < void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
4686 > void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2,
4687 >                                  flagPair pairFlags, double scaleFactor){
4688  
4689    bool sameObjects = false;
4690    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
# Line 3245 | Line 4697 | void OSUAnalysis::fill2DHistogram(TH2* h
4697  
4698        pairCounter++;
4699  
4700 <      //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;
3251 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
4700 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter).second) continue;
4701  
4702 +      string stringValue = "";
4703        string currentString = parameters.inputVariables.at(0);
4704        string inputVariable = "";
4705        string function = "";
4706 <      if(currentString.find("(")==std::string::npos){
4706 >      if(currentString.find("(")==string::npos){
4707          inputVariable = currentString;// variable to cut on
4708        }
4709        else{
# Line 3261 | Line 4711 | void OSUAnalysis::fill2DHistogram(TH2* h
4711          inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4712          inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4713        }
4714 <      double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
4714 >      double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4715  
4716        currentString = parameters.inputVariables.at(1);
4717        inputVariable = "";
4718        function = "";
4719 <      if(currentString.find("(")==std::string::npos){
4719 >      if(currentString.find("(")==string::npos){
4720          inputVariable = currentString;// variable to cut on
4721        }
4722        else{
# Line 3274 | Line 4724 | void OSUAnalysis::fill2DHistogram(TH2* h
4724          inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
4725          inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
4726        }
4727 <      double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
4727 >      double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function, stringValue);
4728  
4729  
4730        histo->Fill(valueX,valueY,scaleFactor);
# Line 3295 | Line 4745 | int OSUAnalysis::getGenMatchedParticleIn
4745  
4746      double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
4747      if(currentDeltaR > 0.05) continue;
4748 <    //     cout << std::setprecision(3) << std::setw(20)
4748 >    //     clog << setprecision(3) << setw(20)
4749      //          << "\tcurrentParticle:  eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
4750 <    //          << std::setw(20)
4750 >    //          << setw(20)
4751      //          << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
4752 <    //          << std::setw(20)
4752 >    //          << setw(20)
4753      //          << "\tdeltaR = " << currentDeltaR
4754 <    //          << std::setprecision(1)
4755 <    //          << std::setw(20)
4754 >    //          << setprecision(1)
4755 >    //          << setw(20)
4756      //          << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
4757 <    //          << std::setw(20)
4757 >    //          << setw(20)
4758      //          << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
4759 <    //          << std::setw(20)
4759 >    //          << setw(20)
4760      //          << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
4761      if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
4762        bestMatchIndex = mcparticle - mcparticles->begin();
# Line 3314 | Line 4764 | int OSUAnalysis::getGenMatchedParticleIn
4764      }
4765  
4766    }
4767 <  //   if(bestMatchDeltaR != 999)  cout << "bestMatch:  deltaR = " << bestMatchDeltaR << "   id = " << mcparticles->at(bestMatchIndex).id << "   motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
4768 <  //   else cout << "no match found..." << endl;
4767 >  //   if(bestMatchDeltaR != 999)  clog << "bestMatch:  deltaR = " << bestMatchDeltaR << "   id = " << mcparticles->at(bestMatchIndex).id << "   motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
4768 >  //   else clog << "no match found..." << endl;
4769    return bestMatchIndex;
4770  
4771   }
# Line 3369 | Line 4819 | int OSUAnalysis::findTauMotherIndex(cons
4819   // 20        strange baryon
4820   // 21        charm baryon
4821   // 22        bottom baryon
4822 < // 23        other
4822 > // 23        QCD string
4823 > // 24        other
4824  
4825   int OSUAnalysis::getPdgIdBinValue(int pdgId){
4826  
# Line 3393 | Line 4844 | int OSUAnalysis::getPdgIdBinValue(int pd
4844    else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
4845    else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
4846    else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
4847 +  else if(absPdgId == 92  ) binValue = 23;
4848  
4849 <  else binValue = 23;
4849 >  else binValue = 24;
4850  
4851    return binValue;
4852  
# Line 3404 | Line 4856 | const BNprimaryvertex *
4856   OSUAnalysis::chosenVertex ()
4857   {
4858    const BNprimaryvertex *chosenVertex = 0;
4859 <  if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
4860 <    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back().size() ? cumulativeFlags.at("primaryvertexs").back() :
4861 <      cumulativeFlags.at("primaryvertexs").at(cumulativeFlags.at("primaryvertexs").size() - 2);
4859 >  if(cumulativeFlags.find ("primaryvertexs") != cumulativeFlags.end ()){
4860 >    flagPair vertexFlags;
4861 >    for (int i = cumulativeFlags.at("primaryvertexs").size() - 1; i >= 0; i--){
4862 >      if (cumulativeFlags.at("primaryvertexs").at(i).size()){
4863 >        vertexFlags = cumulativeFlags.at("primaryvertexs").at(i);
4864 >        break;
4865 >      }
4866 >    }
4867      for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
4868 <      if(!vertexFlags.at(vertexIndex)) continue;
4868 >      if(!vertexFlags.at(vertexIndex).first) continue;
4869        chosenVertex = & primaryvertexs->at(vertexIndex);
4870        break;
4871      }
4872    }
4873 <  else {
4874 <    chosenVertex = &primaryvertexs->at(0);
3418 <  }
4873 >  else if (find (objectsToGet.begin (), objectsToGet.end (), "primaryvertexs") != objectsToGet.end ())
4874 >    chosenVertex = & primaryvertexs->at (0);
4875  
4876    return chosenVertex;
4877   }
# Line 3424 | Line 4880 | const BNmet *
4880   OSUAnalysis::chosenMET ()
4881   {
4882    const BNmet *chosenMET = 0;
4883 <  if(std::find(objectsToCut.begin(), objectsToCut.end(), "mets") != objectsToCut.end()) {
4884 <    vector<bool> metFlags = cumulativeFlags.at("mets").back().size() ? cumulativeFlags.at("mets").back() :
4885 <      cumulativeFlags.at("mets").at(cumulativeFlags.at("mets").size() - 2);
4883 >  if(cumulativeFlags.find ("mets") != cumulativeFlags.end ()){
4884 >    flagPair metFlags;
4885 >    for (int i = cumulativeFlags.at("mets").size() - 1; i >= 0; i--){
4886 >      if (cumulativeFlags.at("mets").at(i).size()){
4887 >        metFlags = cumulativeFlags.at("mets").at(i);
4888 >        break;
4889 >      }
4890 >    }
4891      for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
4892 <      if(!metFlags.at(metIndex)) continue;
4892 >      if(!metFlags.at(metIndex).first) continue;
4893        chosenMET = & mets->at(metIndex);
4894        break;
4895      }
4896    }
4897 <  else {
4898 <    chosenMET = &mets->at(0);
3438 <  }
4897 >  else if (find (objectsToGet.begin (), objectsToGet.end (), "mets") != objectsToGet.end ())
4898 >    chosenMET = & mets->at (0);
4899  
4900    return chosenMET;
4901   }
# Line 3444 | Line 4904 | const BNelectron *
4904   OSUAnalysis::chosenElectron ()
4905   {
4906    const BNelectron *chosenElectron = 0;
4907 <  if(std::find(objectsToCut.begin(), objectsToCut.end(), "electrons") != objectsToCut.end()) {
4908 <    vector<bool> electronFlags = cumulativeFlags.at("electrons").back().size() ? cumulativeFlags.at("electrons").back() :
4909 <      cumulativeFlags.at("electrons").at(cumulativeFlags.at("electrons").size() - 2);
4907 >  if(cumulativeFlags.find ("electrons") != cumulativeFlags.end ()){
4908 >    flagPair electronFlags;
4909 >    for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
4910 >      if (cumulativeFlags.at("electrons").at(i).size()){
4911 >        electronFlags = cumulativeFlags.at("electrons").at(i);
4912 >        break;
4913 >      }
4914 >    }
4915      for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
4916 <      if(!electronFlags.at(electronIndex)) continue;
4916 >      if(!electronFlags.at(electronIndex).first) continue;
4917        chosenElectron = & electrons->at(electronIndex);
4918        break;
4919      }
4920    }
4921 <  else {
4922 <    chosenElectron = &electrons->at(0);
3458 <  }
4921 >  else if (find (objectsToGet.begin (), objectsToGet.end (), "electrons") != objectsToGet.end ())
4922 >    chosenElectron = & electrons->at (0);
4923  
4924    return chosenElectron;
4925   }
4926  
4927 +
4928   const BNmuon *
4929   OSUAnalysis::chosenMuon ()
4930   {
4931    const BNmuon *chosenMuon = 0;
4932 <  if(std::find(objectsToCut.begin(), objectsToCut.end(), "muons") != objectsToCut.end()) {
4933 <    vector<bool> muonFlags = cumulativeFlags.at("muons").back().size() ? cumulativeFlags.at("muons").back() :
4934 <      cumulativeFlags.at("muons").at(cumulativeFlags.at("muons").size() - 2);
4932 >  if(cumulativeFlags.find ("muons") != cumulativeFlags.end ()){
4933 >    flagPair muonFlags;
4934 >    for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
4935 >      if (cumulativeFlags.at("muons").at(i).size()){
4936 >        muonFlags = cumulativeFlags.at("muons").at(i);
4937 >        break;
4938 >      }
4939 >    }
4940      for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
4941 <      if(!muonFlags.at(muonIndex)) continue;
4941 >      if(!muonFlags.at(muonIndex).first) continue;
4942        chosenMuon = & muons->at(muonIndex);
4943        break;
4944      }
4945    }
4946 <  else {
4947 <    chosenMuon = &muons->at(0);
3478 <  }
4946 >  else if (find (objectsToGet.begin (), objectsToGet.end (), "muons") != objectsToGet.end ())
4947 >    chosenMuon = & muons->at (0);
4948  
4949    return chosenMuon;
4950   }
4951  
4952 + double
4953 + OSUAnalysis::chosenHT ()
4954 + {
4955 +  double chosenHT = 0.0;
4956 +  if(cumulativeFlags.find ("jets") != cumulativeFlags.end ()){
4957 +    flagPair jetFlags;
4958 +    for (int i = cumulativeFlags.at("jets").size() - 1; i >= 0; i--){
4959 +      if (cumulativeFlags.at("jets").at(i).size()){
4960 +        jetFlags = cumulativeFlags.at("jets").at(i);
4961 +        break;
4962 +      }
4963 +    }
4964 +    for (uint jetIndex = 0; jetIndex != jetFlags.size(); jetIndex++){
4965 +      if(!jetFlags.at(jetIndex).first) continue;
4966 +      chosenHT += jets->at(jetIndex).pt;
4967 +    }
4968 +  }
4969 +
4970 +  return chosenHT;
4971 + }
4972 +
4973 + pair<const BNmuon *, const BNmuon *>
4974 + OSUAnalysis::leadMuonPair ()
4975 + {
4976 +  pair<const BNmuon *, const BNmuon *> leadMuonPair;
4977 +  leadMuonPair.first = leadMuonPair.second = 0;
4978 +
4979 +  if(cumulativeFlags.find ("muons") != cumulativeFlags.end ()){
4980 +    flagPair muonFlags;
4981 +    for (int i = cumulativeFlags.at("muons").size() - 1; i >= 0; i--){
4982 +      if (cumulativeFlags.at("muons").at(i).size()){
4983 +        muonFlags = cumulativeFlags.at("muons").at(i);
4984 +        break;
4985 +      }
4986 +    }
4987 +    for (uint muonIndex0 = 0; muonIndex0 != muonFlags.size(); muonIndex0++){
4988 +      if(!muonFlags.at(muonIndex0).first) continue;
4989 +      for (uint muonIndex1 = muonIndex0 + 1; muonIndex1 < muonFlags.size(); muonIndex1++){
4990 +        if(!muonFlags.at(muonIndex1).first) continue;
4991 +        const BNmuon *mu0 = & muons->at(muonIndex0);
4992 +        const BNmuon *mu1 = & muons->at(muonIndex1);
4993 +        if(leadMuonPair.first == 0 || leadMuonPair.second == 0){
4994 +          leadMuonPair.first = mu0;
4995 +          leadMuonPair.second = mu1;
4996 +        }
4997 +        else{
4998 +          TVector2 newPt0 (mu0->px, mu0->py),
4999 +                   newPt1 (mu1->px, mu1->py),
5000 +                   oldPt0 (leadMuonPair.first->px, leadMuonPair.first->py),
5001 +                   oldPt1 (leadMuonPair.second->px, leadMuonPair.second->py);
5002 +          if(newPt0.Mod () + newPt1.Mod () > oldPt0.Mod() + oldPt1.Mod ())
5003 +            {
5004 +              leadMuonPair.first = mu0;
5005 +              leadMuonPair.second = mu1;
5006 +            }
5007 +        }
5008 +      }
5009 +    }
5010 +  }
5011 +
5012 +  return leadMuonPair;
5013 + }
5014 +
5015 + pair<const BNelectron *, const BNelectron *>
5016 + OSUAnalysis::leadElectronPair ()
5017 + {
5018 +  pair<const BNelectron *, const BNelectron *> leadElectronPair;
5019 +  leadElectronPair.first = leadElectronPair.second = 0;
5020 +  if(cumulativeFlags.find ("electrons") != cumulativeFlags.end ()){
5021 +    flagPair electronFlags;
5022 +    for (int i = cumulativeFlags.at("electrons").size() - 1; i >= 0; i--){
5023 +      if (cumulativeFlags.at("electrons").at(i).size()){
5024 +        electronFlags = cumulativeFlags.at("electrons").at(i);
5025 +        break;
5026 +      }
5027 +    }
5028 +    for (uint electronIndex0 = 0; electronIndex0 != electronFlags.size(); electronIndex0++){
5029 +      if(!electronFlags.at(electronIndex0).first) continue;
5030 +      for (uint electronIndex1 = electronIndex0 + 1; electronIndex1 < electronFlags.size(); electronIndex1++){
5031 +        if(!electronFlags.at(electronIndex1).first) continue;
5032 +        const BNelectron *el0 = & electrons->at(electronIndex0);
5033 +        const BNelectron *el1 = & electrons->at(electronIndex1);
5034 +        if(leadElectronPair.first == 0 || leadElectronPair.second == 0){
5035 +          leadElectronPair.first = el0;
5036 +          leadElectronPair.second = el1;
5037 +        }
5038 +        else{
5039 +          TVector2 newPt0 (el0->px, el0->py),
5040 +                   newPt1 (el1->px, el1->py),
5041 +                   oldPt0 (leadElectronPair.first->px, leadElectronPair.first->py),
5042 +                   oldPt1 (leadElectronPair.second->px, leadElectronPair.second->py);
5043 +          if(newPt0.Mod () + newPt1.Mod () > oldPt0.Mod() + oldPt1.Mod ())
5044 +            {
5045 +              leadElectronPair.first = el0;
5046 +              leadElectronPair.second = el1;
5047 +            }
5048 +        }
5049 +      }
5050 +    }
5051 +  }
5052 +
5053 +  return leadElectronPair;
5054 + }
5055  
5056   DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines