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.3 by lantonel, Mon Jan 28 13:30:45 2013 UTC vs.
Revision 1.34 by wulsin, Wed Mar 20 20:18:43 2013 UTC

# Line 16 | 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")),
27 +  channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
28 +  histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
29 +  plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
30 +  doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting"))
31 + {
32  
33 +  TH1::SetDefaultSumw2 ();
34  
35 <  channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels"))
35 >  cout << "Debug:  about to start doing pile-up reweighting." << endl;  
36  
37 < {
38 <  TH1::SetDefaultSumw2 ();
37 >  //create pile-up reweighting object, if necessary
38 >  if(datasetType_ != "data") {
39 >    if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
40 >    cout << "Debug PUReweighting 0." << endl;  
41 >    muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
42 >    cout << "Debug PUReweighting 1." << endl;  
43 >    electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
44 >    cout << "Debug PUReweighting 2." << endl;  
45 >  }
46 >
47 >  cout << "Finished pile-up reweighting." << endl;  
48  
49    // Construct Cutflow Objects. These store the results of cut decisions and
50    // handle filling cut flow histograms.
51    masterCutFlow_ = new CutFlow (fs_);
52    std::vector<TFileDirectory> directories;
53  
54 +  //always get vertex collection so we can assign the primary vertex in the event
55 +  objectsToGet.push_back("primaryvertexs");
56 +
57 +  //always make the plot of number of primary verticex (to check pile-up reweighting)
58 +  objectsToPlot.push_back("primaryvertexs");
59 +
60 +  //always get the MC particles to do GEN-matching
61 +  objectsToGet.push_back("mcparticles");
62 +
63 +  //always get the event collection to do pile-up reweighting
64 +  objectsToGet.push_back("events");
65 +
66 +  //parse the histogram definitions
67 +  for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
68 +
69 +    string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
70 +    if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
71 +    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
72 +      objectsToGet.push_back(tempInputCollection);
73 +      objectsToPlot.push_back(tempInputCollection);
74 +      objectsToCut.push_back(tempInputCollection);
75 +    }
76 +    else{//pair of objects, need to add them both to the things to objectsToGet
77 +      int dashIndex = tempInputCollection.find("-");
78 +      int spaceIndex = tempInputCollection.find(" ");
79 +      int secondWordLength = spaceIndex - dashIndex;
80 +      objectsToGet.push_back(tempInputCollection);
81 +      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
82 +      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
83 +      objectsToPlot.push_back(tempInputCollection);
84 +      objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
85 +      objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
86 +      objectsToCut.push_back(tempInputCollection);
87 +      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
88 +      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
89 +
90 +      }
91 +
92 +    vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
93 +
94 +    for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
95 +
96 +      histogram tempHistogram;
97 +      tempHistogram.inputCollection = tempInputCollection;
98 +      tempHistogram.name = histogramList_.at(currentHistogram).getParameter<string>("name");
99 +      tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
100 +      tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
101 +      tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
102 +
103 +      histograms.push_back(tempHistogram);
104 +
105 +    }
106 +  }
107 +  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
108 +  sort( objectsToPlot.begin(), objectsToPlot.end() );
109 +  objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
110 +
111 +
112 +
113 +  //add histograms with the gen-matched id, mother id, and grandmother id
114 +  for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
115 +
116 +    string currentObject = objectsToPlot.at(currentObjectIndex);
117 +    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
118 +
119 +    histogram tempIdHisto;
120 +    histogram tempMomIdHisto;
121 +    histogram tempGmaIdHisto;
122 +
123 +    tempIdHisto.inputCollection = currentObject;
124 +    tempMomIdHisto.inputCollection = currentObject;
125 +    tempGmaIdHisto.inputCollection = currentObject;
126 +
127 +    currentObject = currentObject.substr(0, currentObject.size()-1);
128 +    tempIdHisto.name = currentObject+"GenMatchId";
129 +    tempMomIdHisto.name = currentObject+"GenMatchMotherId";
130 +    tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
131 +
132 +    currentObject.at(0) = toupper(currentObject.at(0));
133 +    tempIdHisto.title = currentObject+" Gen-matched Particle";
134 +    tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
135 +    tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
136 +
137 +    int maxNum = 24;
138 +    vector<double> binVector;
139 +    binVector.push_back(maxNum);
140 +    binVector.push_back(0);
141 +    binVector.push_back(maxNum);
142 +
143 +    tempIdHisto.bins = binVector;
144 +    tempIdHisto.inputVariables.push_back("genMatchedId");
145 +    tempMomIdHisto.bins = binVector;
146 +    tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
147 +    tempGmaIdHisto.bins = binVector;
148 +    tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
149 +
150 +    histograms.push_back(tempIdHisto);
151 +    histograms.push_back(tempMomIdHisto);
152 +    histograms.push_back(tempGmaIdHisto);
153 +
154 +  }
155 +
156 +
157  
158    channel tempChannel;
159    //loop over all channels (event selections)
160    for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
161 <    
161 >
162      //get name of channel
163      string channelName  (channels_.at(currentChannel).getParameter<string>("name"));
164      tempChannel.name = channelName;
# Line 45 | Line 171 | OSUAnalysis::OSUAnalysis (const edm::Par
171      if(channels_.at(currentChannel).exists("triggers")){
172        triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
173        for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
174 <        tempChannel.triggers.push_back(triggerNames.at(trigger));
175 <      allNecessaryObjects.push_back("triggers");
174 >        tempChannel.triggers.push_back(triggerNames.at(trigger));
175 >      objectsToGet.push_back("triggers");
176      }
177  
178  
179 +
180 +
181      //create cutFlow for this channel
182      cutFlows_.push_back (new CutFlow (fs_, channelName));
183  
184      //book a directory in the output file with the name of the channel
185      TFileDirectory subDir = fs_->mkdir( channelName );
186      directories.push_back(subDir);
187 <    std::map<std::string, TH1D*> histoMap;
188 <    oneDHists_.push_back(histoMap);
187 >
188 >    std::map<std::string, TH1D*> oneDhistoMap;
189 >    oneDHists_.push_back(oneDhistoMap);
190 >    std::map<std::string, TH2D*> twoDhistoMap;
191 >    twoDHists_.push_back(twoDhistoMap);
192 >
193 >
194 >    //book all histograms included in the configuration
195 >    for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
196 >      histogram currentHistogram = histograms.at(currentHistogramIndex);
197 >      int numBinsElements = currentHistogram.bins.size();
198 >      int numInputVariables = currentHistogram.inputVariables.size();
199 >
200 >      if(numBinsElements != 3 && numBinsElements !=6) {
201 >        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
202 >        exit(0);
203 >      }
204 >      else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
205 >        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
206 >        exit(0);
207 >      }
208 >      else if(numBinsElements == 3){
209 >        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));
210 >      }
211 >      else if(numBinsElements == 6){
212 >        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));
213 >
214 >      }
215 >
216 >
217 >      if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
218 >
219 > // bin      particle type
220 > // ---      -------------
221 > //  0        unmatched
222 > //  1        u
223 > //  2        d
224 > //  3        s
225 > //  4        c
226 > //  5        b
227 > //  6        t
228 > //  7        e
229 > //  8        mu
230 > //  9        tau
231 > // 10        nu
232 > // 11        g
233 > // 12        gamma
234 > // 13        Z
235 > // 14        W
236 > // 15        light meson
237 > // 16        K meson
238 > // 17        D meson
239 > // 18        B meson
240 > // 19        light baryon
241 > // 20        strange baryon
242 > // 21        charm baryon
243 > // 22        bottom baryon
244 > // 23        other
245 >
246 >      vector<TString> labelArray;
247 >      labelArray.push_back("unmatched");
248 >      labelArray.push_back("u");
249 >      labelArray.push_back("d");
250 >      labelArray.push_back("s");
251 >      labelArray.push_back("c");
252 >      labelArray.push_back("b");
253 >      labelArray.push_back("t");
254 >      labelArray.push_back("e");
255 >      labelArray.push_back("#mu");
256 >      labelArray.push_back("#tau");
257 >      labelArray.push_back("#nu");
258 >      labelArray.push_back("g");
259 >      labelArray.push_back("#gamma");
260 >      labelArray.push_back("Z");
261 >      labelArray.push_back("W");
262 >      labelArray.push_back("light meson");
263 >      labelArray.push_back("K meson");
264 >      labelArray.push_back("D meson");
265 >      labelArray.push_back("B meson");
266 >      labelArray.push_back("light baryon");
267 >      labelArray.push_back("strange baryon");
268 >      labelArray.push_back("charm baryon");
269 >      labelArray.push_back("bottom baryon");
270 >      labelArray.push_back("other");
271 >
272 >      for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
273 >        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
274 >      }
275 >    }
276 >
277 >    //book a histogram for the number of each object type to be plotted
278 >
279 >    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
280 >      string currentObject = objectsToPlot.at(currentObjectIndex);
281 >      int maxNum = 10;
282 >      if(currentObject == "mcparticles") maxNum = 50;
283 >      else if(currentObject == "primaryvertexs") maxNum = 50;
284 >      else if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
285 >      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
286 >      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
287 >
288 >      currentObject.at(0) = toupper(currentObject.at(0));
289 >      string histoName = "num" + currentObject;
290 >
291 >      if(histoName == "numPrimaryvertexs"){
292 >        string newHistoName = histoName + "BeforePileupCorrection";
293 >        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
294 >        newHistoName = histoName + "AfterPileupCorrection";
295 >        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
296 >      }
297 >      else
298 >        oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
299 >    }
300  
301  
63    //muon histograms
64    allNecessaryObjects.push_back("muons");
65    oneDHists_.at(currentChannel)["muonPt"]  = directories.at(currentChannel).make<TH1D> ("muonPt",channelLabel+" channel: Muon Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
66    oneDHists_.at(currentChannel)["muonEta"] = directories.at(currentChannel).make<TH1D> ("muonEta",channelLabel+" channel: Muon Eta; #eta", 100, -5, 5);
67    oneDHists_.at(currentChannel)["muonD0"] = directories.at(currentChannel).make<TH1D> ("muonD0",channelLabel+" channel: Muon d_{0}; d_{0} [cm] ", 1000, -1, 1);
68    oneDHists_.at(currentChannel)["muonAbsD0"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0",channelLabel+" channel: Muon d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
69    oneDHists_.at(currentChannel)["muonD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonD0Sig",channelLabel+" channel: Muon d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
70    oneDHists_.at(currentChannel)["muonAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0Sig",channelLabel+" channel: Muon d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
71    oneDHists_.at(currentChannel)["muonIso"] = directories.at(currentChannel).make<TH1D> ("muonIso",channelLabel+" channel: Muon Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
72    oneDHists_.at(currentChannel)["numMuons"] = directories.at(currentChannel).make<TH1D> ("numMuons",channelLabel+" channel: Number of Selected Muons; # muons", 10, 0, 10);
73
74
75    //electron histograms
76    allNecessaryObjects.push_back("electrons");
77    oneDHists_.at(currentChannel)["electronPt"]  = directories.at(currentChannel).make<TH1D> ("electronPt",channelLabel+" channel: Electron Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
78    oneDHists_.at(currentChannel)["electronEta"] = directories.at(currentChannel).make<TH1D> ("electronEta",channelLabel+" channel: Electron Eta; #eta", 50, -5, 5);
79    oneDHists_.at(currentChannel)["electronD0"] = directories.at(currentChannel).make<TH1D> ("electronD0",channelLabel+" channel: Electron d_{0}; d_{0} [cm] ", 1000, -1, 1);
80    oneDHists_.at(currentChannel)["electronAbsD0"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0",channelLabel+" channel: Electron d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
81    oneDHists_.at(currentChannel)["electronD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronD0Sig",channelLabel+" channel: Electron d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
82    oneDHists_.at(currentChannel)["electronAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0Sig",channelLabel+" channel: Electron d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
83    oneDHists_.at(currentChannel)["electronIso"] = directories.at(currentChannel).make<TH1D> ("electronIso",channelLabel+" channel: Electron Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
84    oneDHists_.at(currentChannel)["electronFbrem"] = directories.at(currentChannel).make<TH1D> ("electronFbrem",channelLabel+" channel: Electron Brem. Energy Fraction; fbrem", 1000, 0, 2);
85    oneDHists_.at(currentChannel)["numElectrons"] = directories.at(currentChannel).make<TH1D> ("numElectrons",channelLabel+" channel: Number of Selected Electrons; # electrons", 10, 0, 10);
302  
303  
304  
305      //get list of cuts for this channel
306      vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
307  
92    vector<string> tempInputCollections;
93
94
308      //loop over and parse all cuts
309      for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
97
310        cut tempCut;
311 <     //store input collection for cut
312 <      string inputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
313 <      tempCut.inputCollection = inputCollection;
314 <      
315 <      tempInputCollections.push_back(inputCollection);
316 <      allNecessaryObjects.push_back(inputCollection);
311 >      //store input collection for cut
312 >      string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
313 >      tempCut.inputCollection = tempInputCollection;
314 >      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
315 >        objectsToGet.push_back(tempInputCollection);
316 >        objectsToCut.push_back(tempInputCollection);
317 >      }
318 >      else{//pair of objects, need to add them both to the things to objectsToGet
319 >        int dashIndex = tempInputCollection.find("-");
320 >        int spaceIndex = tempInputCollection.find(" ");
321 >        int secondWordLength = spaceIndex - dashIndex;
322 >        objectsToGet.push_back(tempInputCollection);
323 >        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
324 >        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
325 >        objectsToCut.push_back(tempInputCollection);
326 >        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
327 >        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
328 >
329 >      }
330 >
331 >
332  
333        //split cut string into parts and store them
334        string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
335        std::vector<string> cutStringVector = splitString(cutString);
336 <      tempCut.variable = cutStringVector.at(0);// variable to cut on
337 <      tempCut.comparativeOperator = cutStringVector.at(1);// comparison to make
338 <      tempCut.cutValue = atof(cutStringVector.at(2).c_str());// threshold value to pass cut
336 >      if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
337 >        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
338 >        exit(0);
339 >      }
340 >      tempCut.numSubcuts = (cutStringVector.size()+1)/4;
341 >      for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
342 >        int indexOffset = 4 * subcutIndex;
343 >        string currentVariableString = cutStringVector.at(indexOffset);
344 >        if(currentVariableString.find("(")==std::string::npos){
345 >          tempCut.functions.push_back("");//no function was specified
346 >          tempCut.variables.push_back(currentVariableString);// variable to cut on
347 >        }
348 >        else{
349 >          tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
350 >          string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
351 >          tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
352 >        }
353 >        tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
354 >        tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
355 >        if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
356 >      }
357  
358        //get number of objects required to pass cut for event to pass
359        string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
360        std::vector<string> numberRequiredVector = splitString(numberRequiredString);
361 +      if(numberRequiredVector.size()!=2){
362 +        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
363 +        exit(0);
364 +      }
365  
366 <      // determine number required if comparison contains "="
118 <      int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
119 <      if(numberRequiredVector.at(0) == ">") numberRequiredInt++;
120 <      else if(numberRequiredVector.at(0) == "<") numberRequiredInt--;
121 <
366 >      int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
367        tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
368        tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
369 <      
370 <      if(cuts_.at(currentCut).exists("function")){
126 <        tempCut.function = cuts_.at(currentCut).getParameter<string> ("function");
127 <      }
128 <      else tempCut.function = "";
369 >
370 >
371        string tempCutName;
372        if(cuts_.at(currentCut).exists("alias")){
373 <        tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
373 >        tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
374        }
375        else{
376 <        //construct string for cutflow table
377 <        bool plural = numberRequiredInt != 1;
378 <        string collectionString = plural ? inputCollection : inputCollection.substr(0, inputCollection.size()-1);
379 <        string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
380 <        tempCutName = cutName;
376 >        //construct string for cutflow table
377 >        bool plural = numberRequiredInt != 1;
378 >        string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
379 >        string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
380 >        tempCutName = cutName;
381        }
382        tempCut.name = tempCutName;
383  
# Line 143 | Line 385 | OSUAnalysis::OSUAnalysis (const edm::Par
385        tempChannel.cuts.push_back(tempCut);
386  
387  
146    }//end loop over cuts
388  
389 <    //make unique vector of all objects that this channel cuts on (so we know to set flags for those)
149 <    sort( tempInputCollections.begin(), tempInputCollections.end() );
150 <    tempInputCollections.erase( unique( tempInputCollections.begin(), tempInputCollections.end() ), tempInputCollections.end() );
151 <    tempChannel.inputCollections = tempInputCollections;    
389 >    }//end loop over cuts
390  
391      channels.push_back(tempChannel);
392      tempChannel.cuts.clear();
393  
394    }//end loop over channels
395  
396 <  //make unique vector of objects we need to cut on (so we make sure to get them from the event)
397 <  sort( allNecessaryObjects.begin(), allNecessaryObjects.end() );
398 <  allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() );
396 >
397 >  //make unique vector of objects we need to get from the event
398 >  sort( objectsToGet.begin(), objectsToGet.end() );
399 >  objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
400 >  //make unique vector of objects we need to cut on
401 >  sort( objectsToCut.begin(), objectsToCut.end() );
402 >  objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
403  
404  
405   }
# Line 175 | Line 417 | void
417   OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
418   {
419  
178
420    // Retrieve necessary collections from the event.
421 <  edm::Handle<BNtriggerCollection> triggers;
422 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "triggers") != allNecessaryObjects.end())
421 >
422 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
423      event.getByLabel (triggers_, triggers);
424  
425 <  edm::Handle<BNjetCollection> jets;
185 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end())  
425 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
426      event.getByLabel (jets_, jets);
427  
428 <  edm::Handle<BNmuonCollection> muons;
189 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end())  
428 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
429      event.getByLabel (muons_, muons);
430  
431 <  edm::Handle<BNelectronCollection> electrons;
193 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end())  
431 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
432      event.getByLabel (electrons_, electrons);
433  
434 <  edm::Handle<BNeventCollection> events;
197 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end())  
434 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
435      event.getByLabel (events_, events);
436  
437 <  edm::Handle<BNtauCollection> taus;
201 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end())  
437 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
438      event.getByLabel (taus_, taus);
439  
440 <  edm::Handle<BNmetCollection> mets;
205 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end())  
440 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
441      event.getByLabel (mets_, mets);
442  
443 <  edm::Handle<BNtrackCollection> tracks;
209 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end())  
443 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
444      event.getByLabel (tracks_, tracks);
445  
446 <  edm::Handle<BNgenjetCollection> genjets;
213 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end())  
446 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
447      event.getByLabel (genjets_, genjets);
448  
449 <  edm::Handle<BNmcparticleCollection> mcparticles;
217 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end())  
449 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
450      event.getByLabel (mcparticles_, mcparticles);
451  
452 <  edm::Handle<BNprimaryvertexCollection> primaryvertexs;
221 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end())  
452 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
453      event.getByLabel (primaryvertexs_, primaryvertexs);
454  
455 <  edm::Handle<BNbxlumiCollection> bxlumis;
225 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end())  
455 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
456      event.getByLabel (bxlumis_, bxlumis);
457  
458 <  edm::Handle<BNphotonCollection> photons;
229 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end())  
458 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
459      event.getByLabel (photons_, photons);
460  
461 <  edm::Handle<BNsuperclusterCollection> superclusters;
233 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end())  
461 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
462      event.getByLabel (superclusters_, superclusters);
463  
464 +  //get pile-up event weight
465 +  double scaleFactor = 1.00;
466 +  if(doPileupReweighting_ && datasetType_ != "data")
467 +    scaleFactor = puWeight_->at (events->at (0).numTruePV);
468 +
469  
470    //loop over all channels
471  
472    for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
473      channel currentChannel = channels.at(currentChannelIndex);
241    
242    flagMap individualFlags;  
243    flagMap cumulativeFlags;  
244    counterMap passingCounter;
474  
475 +    flagMap individualFlags;
476 +    counterMap passingCounter;
477 +    cumulativeFlags.clear ();
478 +
479 +    bool triggerDecision = true;
480      if(currentChannel.triggers.size() != 0){  //triggers specified
481 <      bool triggerDecision = evaluateTriggers(currentChannel.triggers,triggers.product());
481 >      triggerDecision = evaluateTriggers(currentChannel.triggers,triggers.product());
482        cutFlows_.at(currentChannelIndex)->at ("trigger") = triggerDecision;
483      }
484  
485      //loop over all cuts
486 +
487 +
488 +
489      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
490        cut currentCut = currentChannel.cuts.at(currentCutIndex);
491  
492 <      for(uint currentObjectIndex = 0; currentObjectIndex != allNecessaryObjects.size(); currentObjectIndex++){
492 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
493 >        string currentObject = objectsToCut.at(currentObjectIndex);
494 >
495 >        //initialize maps to get ready to set cuts
496 >        individualFlags[currentObject].push_back (vector<bool> ());
497 >        cumulativeFlags[currentObject].push_back (vector<bool> ());
498 >
499 >      }
500 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
501 >        string currentObject = objectsToCut.at(currentObjectIndex);
502 >
503 >        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
504 >
505  
506 <        string currentObject = allNecessaryObjects.at(currentObjectIndex);
507 <        individualFlags[currentObject].push_back (vector<bool> ());
508 <        cumulativeFlags[currentObject].push_back (vector<bool> ());
509 <
510 <        if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
511 <        else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
512 <        else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
513 <        else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
514 <        else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
515 <        else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
516 <        else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
517 <        else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
518 <        else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
519 <        else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
520 <        else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
521 <        else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
522 <        else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
506 >        if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
507 >        else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
508 >        else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
509 >        else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
510 >        else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
511 >        else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
512 >        else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
513 >        else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
514 >        else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
515 >        else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
516 >        else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
517 >        else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
518 >        else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
519 >
520 >
521 >        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
522 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
523 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
524 >                                                                   "muon-muon pairs");
525 >        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
526 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
527 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
528 >                                                                           "electron-electron pairs");
529 >        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
530 >                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
531 >                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
532 >                                                                       "electron-muon pairs");
533  
534  
535        }
# Line 278 | Line 537 | OSUAnalysis::analyze (const edm::Event &
537  
538  
539      }//end loop over all cuts
281    
282
540  
541  
542      //use cumulative flags to apply cuts at event level
543  
544      bool eventPassedAllCuts = true;
545  
546 +    //apply trigger (true if none were specified)
547 +    eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
548 +
549 +
550 +
551      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
552  
553        //loop over all objects and count how many passed the cumulative selection up to this point
554        cut currentCut = currentChannel.cuts.at(currentCutIndex);
555        int numberPassing = 0;
556  
557 <      for (uint object = 0; object != cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).size() ; object++)
558 <          if(cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).at(object)) numberPassing++;
559 <
557 >      for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
558 >          if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
559 >      }
560  
561        bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
562        cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
# Line 303 | Line 565 | OSUAnalysis::analyze (const edm::Event &
565  
566      }
567  
568 <    cutFlows_.at(currentChannelIndex)->fillCutFlow();  
569 <  
568 >    cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
569 >
570  
571  
572      if(!eventPassedAllCuts)continue;
573  
574  
575 +    //if(datasetType_ != "data") {
576 +    //  scaleFactor *= muonSFWeight_->at (chosenMuon ()->eta);
577 +    //  scaleFactor *= electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
578 +    //}
579 +
580 +    //filling histograms
581 +    for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
582 +      histogram currentHistogram = histograms.at(histogramIndex);
583 +
584 +      if(currentHistogram.inputVariables.size() == 1){
585 +        TH1D* histo;
586 +        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
587  
314    vector<bool> electronFlags = cumulativeFlags["electrons"].back();
315    int electronCounter = 0;
316
317    for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
318      if(!electronFlags.at(electronIndex)) continue;
588  
589  
590 <      electronCounter++;
591 <      oneDHists_.at(currentChannelIndex)["electronPt"]->Fill (electrons->at(electronIndex).pt);
592 <      oneDHists_.at(currentChannelIndex)["electronEta"]->Fill (electrons->at(electronIndex).eta);
593 <      oneDHists_.at(currentChannelIndex)["electronD0"]->Fill (electrons->at(electronIndex).correctedD0Vertex);
594 <      oneDHists_.at(currentChannelIndex)["electronAbsD0"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex));
595 <      oneDHists_.at(currentChannelIndex)["electronD0Sig"]->Fill (electrons->at(electronIndex).correctedD0Vertex / electrons->at(electronIndex).tkD0err);
596 <      oneDHists_.at(currentChannelIndex)["electronAbsD0Sig"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex) / electrons->at(electronIndex).tkD0err);
597 <      oneDHists_.at(currentChannelIndex)["electronIso"]->Fill (((electrons->at(electronIndex).trackIso + electrons->at(electronIndex).caloIso) / electrons->at(electronIndex).pt));
598 <      oneDHists_.at(currentChannelIndex)["electronFbrem"]->Fill (electrons->at(electronIndex).fbrem);
599 <    }
600 <    oneDHists_.at(currentChannelIndex)["numElectrons"]->Fill (electronCounter);
601 <
602 <    
603 <    vector<bool> muonFlags = cumulativeFlags["muons"].back();
604 <    int muonCounter = 0;
590 >        if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
591 >         else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
592 >        else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
593 >                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
594 >                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
595 >        else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
596 >        else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
597 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
598 >                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
599 >        else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
600 >                                                                                              cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
601 >                                                                                              cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
602 >        else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
603 >        else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
604 >        else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
605 >        else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
606 >        else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
607 >        else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
608 >        else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
609 >        else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
610 >        else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
611 >        else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
612 >      }
613 >      else if(currentHistogram.inputVariables.size() == 2){
614 >        TH2D* histo;
615 >        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
616  
617  
338    for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
339      if(!muonFlags.at(muonIndex)) continue;
340      muonCounter++;
618  
619 <      oneDHists_.at(currentChannelIndex)["muonPt"]->Fill (muons->at(muonIndex).pt);
620 <      oneDHists_.at(currentChannelIndex)["muonEta"]->Fill (muons->at(muonIndex).eta);
621 <      oneDHists_.at(currentChannelIndex)["muonD0"]->Fill (muons->at(muonIndex).correctedD0Vertex);
622 <      oneDHists_.at(currentChannelIndex)["muonAbsD0"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex));
623 <      oneDHists_.at(currentChannelIndex)["muonD0Sig"]->Fill (muons->at(muonIndex).correctedD0Vertex / muons->at(muonIndex).tkD0err);
624 <      oneDHists_.at(currentChannelIndex)["muonAbsD0Sig"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex) / muons->at(muonIndex).tkD0err);
625 <      oneDHists_.at(currentChannelIndex)["muonIso"]->Fill (((muons->at(muonIndex).trackIso + muons->at(muonIndex).caloIso) / muons->at(muonIndex).pt));
619 >        if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
620 >        else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
621 >        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
622 >                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
623 >                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
624 >        else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
625 >        else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
626 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
627 >                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
628 >        else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
629 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
630 >                                                                                               cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
631 >        else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
632 >        else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
633 >        else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
634 >        else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
635 >        else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
636 >        else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
637 >        else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
638 >        else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
639 >        else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
640 >        else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
641 >      }
642      }
643 <    oneDHists_.at(currentChannelIndex)["numMuons"]->Fill (muonCounter);
643 >
644 >    //fills histograms with the sizes of collections
645 >    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
646 >      string currentObject = objectsToPlot.at(currentObjectIndex);
647 >
648 >      string objectToPlot = "";
649 >
650 >      if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
651 >      else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
652 >      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
653 >      else objectToPlot = currentObject;
654 >      string tempCurrentObject = objectToPlot;
655 >      tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
656 >      string histoName = "num" + tempCurrentObject;
657  
658  
659  
660  
661 +      //set position of primary vertex in event, in order to calculate quantities relative to it
662 +      if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
663 +        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
664 +        int numToPlot = 0;
665 +        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
666 +          if(lastCutFlags.at(currentFlag)) numToPlot++;
667 +        }
668 +        if(objectToPlot == "primaryvertexs"){
669 +          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
670 +          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
671 +        }
672 +        else
673 +          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,scaleFactor);
674 +      }
675 +      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),scaleFactor);
676 +      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),scaleFactor);
677 +      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,scaleFactor);
678 +      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),scaleFactor);
679 +      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,scaleFactor);
680 +      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),scaleFactor);
681 +      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),scaleFactor);
682 +      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),scaleFactor);
683 +      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),scaleFactor);
684 +      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),scaleFactor);
685 +      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),scaleFactor);
686 +      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),scaleFactor);
687 +      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),scaleFactor);
688 +      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),scaleFactor);
689 +      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),scaleFactor);
690 +      else if(objectToPlot == "primaryvertexs"){
691 +        oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
692 +        oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
693 +      }
694 +
695 +    }
696 +
697 +
698 +
699  
700    } //end loop over channel
701  
702 <  masterCutFlow_->fillCutFlow();    
702 >  masterCutFlow_->fillCutFlow(scaleFactor);
703 >
704  
705  
706   }
707  
708  
709 < bool
709 > bool
710   OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
711  
712  
# Line 369 | Line 714 | OSUAnalysis::evaluateComparison (double
714    else if(comparison == ">=") return testValue >= cutValue;
715    else if(comparison == "<")  return testValue <  cutValue;
716    else if(comparison == "<=") return testValue <= cutValue;
717 +  else if(comparison == "==") return testValue == cutValue;
718 +  else if(comparison == "=") return testValue == cutValue;
719 +  else if(comparison == "!=") return testValue != cutValue;
720    else {std::cout << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
721  
722   }
# Line 384 | Line 732 | OSUAnalysis::evaluateTriggers (vector<st
732        for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++)
733          {
734            if(trigger->name.find(triggersToTest.at(triggerName))!=std::string::npos){
735 <            triggerDecision = true;
735 >            triggerDecision = true;
736            }
737          }
738      }
# Line 522 | Line 870 | OSUAnalysis::valueLookup (const BNjet* o
870    else if(variable == "puJetId_loose_full") value = object->puJetId_loose_full;
871    else if(variable == "puJetId_loose_simple") value = object->puJetId_loose_simple;
872    else if(variable == "puJetId_loose_cutbased") value = object->puJetId_loose_cutbased;
873 <  
874 <  
873 >
874 >
875    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
876 <  
877 <  if(function == "abs") value = abs(value);
876 >
877 >  value = applyFunction(function, value);
878  
879    return value;
880   }
# Line 550 | Line 898 | OSUAnalysis::valueLookup (const BNmuon*
898    else if(variable == "ecalIso") value = object->ecalIso;
899    else if(variable == "hcalIso") value = object->hcalIso;
900    else if(variable == "caloIso") value = object->caloIso;
901 <  else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
901 >  else if(variable == "trackIsDR03") value = object->trackIsoDR03;
902    else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
903    else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
904    else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
# Line 688 | Line 1036 | OSUAnalysis::valueLookup (const BNmuon*
1036    else if(variable == "numberOfMatchedStations") value = object->numberOfMatchedStations;
1037    else if(variable == "time_ndof") value = object->time_ndof;
1038  
1039 +  //user-defined variables
1040 +  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1041 +  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1042 +  else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1043 +  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1044 +  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1045 +  else if(variable == "metMT") {
1046 +    const BNmet *met = chosenMET ();
1047 +    if (met)
1048 +      {
1049 +        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1050 +                       p2 (met->px, met->py, 0.0, met->pt);
1051 +
1052 +        value = (p1 + p2).Mt ();
1053 +      }
1054 +    else
1055 +      value = -999;
1056 +  }
1057 +
1058 +
1059 +
1060 +  else if(variable == "correctedD0VertexInEBPlus"){
1061 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1062 +    else value = -999;
1063 +  }
1064 +  else if(variable == "correctedD0VertexOutEBPlus"){
1065 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1066 +    else value = -999;
1067 +  }
1068 +  else if(variable == "correctedD0VertexEEPlus"){
1069 +    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex;
1070 +    else value = -999;
1071 +  }
1072 +
1073 +  else if(variable == "correctedD0BeamspotInEBPlus"){
1074 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1075 +    else value = -999;
1076 +  }
1077 +  else if(variable == "correctedD0BeamspotOutEBPlus"){
1078 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1079 +    else value = -999;
1080 +  }
1081 +  else if(variable == "correctedD0BeamspotEEPlus"){
1082 +    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0;
1083 +    else value = -999;
1084 +  }
1085 +
1086 +  else if(variable == "correctedD0VertexInEBMinus"){
1087 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1088 +    else value = -999;
1089 +  }
1090 +  else if(variable == "correctedD0VertexOutEBMinus"){
1091 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1092 +    else value = -999;
1093 +  }
1094 +  else if(variable == "correctedD0VertexEEMinus"){
1095 +    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex;
1096 +    else value = -999;
1097 +  }
1098 +
1099 +  else if(variable == "correctedD0BeamspotInEBMinus"){
1100 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1101 +    else value = -999;
1102 +  }
1103 +  else if(variable == "correctedD0BeamspotOutEBMinus"){
1104 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1105 +    else value = -999;
1106 +  }
1107 +  else if(variable == "correctedD0BeamspotEEMinus"){
1108 +    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0;
1109 +    else value = -999;
1110 +  }
1111 +
1112 +
1113 +  else if(variable == "correctedD0VertexInEBPositiveCharge"){
1114 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1115 +    else value = -999;
1116 +  }
1117 +  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1118 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1119 +    else value = -999;
1120 +  }
1121 +  else if(variable == "correctedD0VertexEEPositiveCharge"){
1122 +    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex;
1123 +    else value = -999;
1124 +  }
1125 +
1126 +  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1127 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1128 +    else value = -999;
1129 +  }
1130 +  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1131 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1132 +    else value = -999;
1133 +  }
1134 +  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1135 +    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0;
1136 +    else value = -999;
1137 +  }
1138 +
1139 +  else if(variable == "correctedD0VertexInEBNegativeCharge"){
1140 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1141 +    else value = -999;
1142 +  }
1143 +  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1144 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1145 +    else value = -999;
1146 +  }
1147 +  else if(variable == "correctedD0VertexEENegativeCharge"){
1148 +    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex;
1149 +    else value = -999;
1150 +  }
1151 +
1152 +  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1153 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1154 +    else value = -999;
1155 +  }
1156 +  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1157 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1158 +    else value = -999;
1159 +  }
1160 +  else if(variable == "correctedD0BeamspotEENegativeCharge"){
1161 +    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0;
1162 +    else value = -999;
1163 +  }
1164 +
1165 +
1166 +  else if(variable == "tightID") {
1167 +    value = object->isGlobalMuon > 0                \
1168 +      && object->isPFMuon > 0                        \
1169 +      && object->normalizedChi2 < 10                \
1170 +      && object->numberOfValidMuonHits > 0        \
1171 +      && object->numberOfMatchedStations > 1        \
1172 +      && fabs(object->correctedD0Vertex) < 0.2        \
1173 +      && fabs(object->correctedDZ) < 0.5        \
1174 +      && object->numberOfValidPixelHits > 0                \
1175 +      && object->numberOfLayersWithMeasurement > 5;
1176 +  }
1177 +  else if(variable == "tightIDdisplaced"){
1178 +    value = object->isGlobalMuon > 0                \
1179 +      && object->normalizedChi2 < 10                \
1180 +      && object->numberOfValidMuonHits > 0        \
1181 +      && object->numberOfMatchedStations > 1        \
1182 +      && object->numberOfValidPixelHits > 0        \
1183 +      && object->numberOfLayersWithMeasurement > 5;
1184 +  }
1185 +
1186 +
1187 +
1188 +  else if(variable == "genMatchedId"){
1189 +    int index = getGenMatchedParticleIndex(object);
1190 +    if(index == -1) value = 0;
1191 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1192 +  }
1193 +  else if(variable == "genMatchedMotherId"){
1194 +    int index = getGenMatchedParticleIndex(object);
1195 +    if(index == -1) value = 0;
1196 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1197 +  }
1198 +  else if(variable == "genMatchedGrandmotherId"){
1199 +    int index = getGenMatchedParticleIndex(object);
1200 +    if(index == -1) value = 0;
1201 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1202 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1203 +      if(motherIndex == -1) value = 0;
1204 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1205 +    }
1206 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1207 +  }
1208 +
1209 +
1210 +
1211    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1212 <  
1212 >
1213 >  value = applyFunction(function, value);
1214  
1215    return value;
1216   }
# Line 856 | Line 1377 | OSUAnalysis::valueLookup (const BNelectr
1377    else if(variable == "eidHyperTight4MC") value = object->eidHyperTight4MC;
1378    else if(variable == "passConvVeto") value = object->passConvVeto;
1379  
1380 +  //user-defined variables
1381 +  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1382 +  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
1383 +  else if(variable == "detIso") value = (object->trackIso) / object->pt;
1384 +  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1385 +  else if(variable == "metMT") {
1386 +    const BNmet *met = chosenMET ();
1387 +    if (met)
1388 +      {
1389 +        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
1390 +                       p2 (met->px, met->py, 0.0, met->pt);
1391 +
1392 +        value = (p1 + p2).Mt ();
1393 +      }
1394 +    else
1395 +      value = -999;
1396 +  }
1397 +
1398 +  else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
1399 +    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
1400 +    else value = -999;
1401 +  }
1402 +  else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){
1403 +    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex;
1404 +    else value = -999;
1405 +  }
1406 +
1407 +  else if(variable == "correctedD0VertexInEBPlus"){
1408 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1409 +    else value = -999;
1410 +  }
1411 +  else if(variable == "correctedD0VertexOutEBPlus"){
1412 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1413 +    else value = -999;
1414 +  }
1415 +  else if(variable == "correctedD0VertexEEPlus"){
1416 +    if(object->isEE && object->eta > 0) value = object->correctedD0Vertex;
1417 +    else value = -999;
1418 +  }
1419 +
1420 +  else if(variable == "correctedD0BeamspotInEBPlus"){
1421 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1422 +    else value = -999;
1423 +  }
1424 +  else if(variable == "correctedD0BeamspotOutEBPlus"){
1425 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1426 +    else value = -999;
1427 +  }
1428 +  else if(variable == "correctedD0BeamspotEEPlus"){
1429 +    if(object->isEE && object->eta > 0) value = object->correctedD0;
1430 +    else value = -999;
1431 +  }
1432 +
1433 +  else if(variable == "correctedD0VertexInEBMinus"){
1434 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1435 +    else value = -999;
1436 +  }
1437 +  else if(variable == "correctedD0VertexOutEBMinus"){
1438 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1439 +    else value = -999;
1440 +  }
1441 +  else if(variable == "correctedD0VertexEEMinus"){
1442 +    if(object->isEE && object->eta < 0) value = object->correctedD0Vertex;
1443 +    else value = -999;
1444 +  }
1445 +
1446 +  else if(variable == "correctedD0BeamspotInEBMinus"){
1447 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1448 +    else value = -999;
1449 +  }
1450 +  else if(variable == "correctedD0BeamspotOutEBMinus"){
1451 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1452 +    else value = -999;
1453 +  }
1454 +  else if(variable == "correctedD0BeamspotEEMinus"){
1455 +    if(object->isEE && object->eta < 0) value = object->correctedD0;
1456 +    else value = -999;
1457 +  }
1458 +
1459 +  else if(variable == "tightID"){
1460 +    if (object->isEB)
1461 +      {
1462 +        value = fabs(object->delEtaIn) < 0.004 \
1463 +          && fabs (object->delPhiIn) < 0.03 \
1464 +          && object->scSigmaIEtaIEta < 0.01 \
1465 +          && object->hadOverEm < 0.12 \
1466 +          && fabs (object->correctedD0Vertex) < 0.02 \
1467 +          && fabs (object->correctedDZ) < 0.1 \
1468 +          && object->absInvEMinusInvPin < 0.05 \
1469 +          && object->passConvVeto;
1470 +      }
1471 +    else
1472 +      {
1473 +        value = fabs(object->delEtaIn) < 0.005 \
1474 +          && fabs (object->delPhiIn) < 0.02 \
1475 +          && object->scSigmaIEtaIEta < 0.03 \
1476 +          && object->hadOverEm < 0.10 \
1477 +          && fabs (object->correctedD0Vertex) < 0.02 \
1478 +          && fabs (object->correctedDZ) < 0.1 \
1479 +          && object->absInvEMinusInvPin < 0.05 \
1480 +          && object->passConvVeto;
1481 +      }
1482 +  }
1483 +
1484 +  else if(variable == "correctedD0VertexInEBPositiveCharge"){
1485 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1486 +    else value = -999;
1487 +  }
1488 +  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1489 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1490 +    else value = -999;
1491 +  }
1492 +  else if(variable == "correctedD0VertexEEPositiveCharge"){
1493 +    if(object->isEE && object->charge > 0) value = object->correctedD0Vertex;
1494 +    else value = -999;
1495 +  }
1496 +
1497 +  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1498 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1499 +    else value = -999;
1500 +  }
1501 +  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1502 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1503 +    else value = -999;
1504 +  }
1505 +  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1506 +    if(object->isEE && object->charge > 0) value = object->correctedD0;
1507 +    else value = -999;
1508 +  }
1509 +
1510 +  else if(variable == "correctedD0VertexInEBNegativeCharge"){
1511 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1512 +    else value = -999;
1513 +  }
1514 +  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1515 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1516 +    else value = -999;
1517 +  }
1518 +  else if(variable == "correctedD0VertexEENegativeCharge"){
1519 +    if(object->isEE && object->charge < 0) value = object->correctedD0Vertex;
1520 +    else value = -999;
1521 +  }
1522 +
1523 +  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1524 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1525 +    else value = -999;
1526 +  }
1527 +  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1528 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1529 +    else value = -999;
1530 +  }
1531 +  else if(variable == "correctedD0BeamspotEENegativeCharge"){
1532 +    if(object->isEE && object->charge < 0) value = object->correctedD0;
1533 +    else value = -999;
1534 +  }
1535 +
1536 +
1537 +  else if(variable == "tightIDdisplaced"){
1538 +    if (object->isEB)
1539 +      {
1540 +        value = fabs(object->delEtaIn) < 0.004 \
1541 +          && fabs (object->delPhiIn) < 0.03 \
1542 +          && object->scSigmaIEtaIEta < 0.01 \
1543 +          && object->hadOverEm < 0.12 \
1544 +          && object->absInvEMinusInvPin < 0.05;
1545 +      }
1546 +    else
1547 +      {
1548 +        value = fabs (object->delEtaIn) < 0.005 \
1549 +          && fabs (object->delPhiIn) < 0.02 \
1550 +          && object->scSigmaIEtaIEta < 0.03 \
1551 +          && object->hadOverEm < 0.10 \
1552 +          && object->absInvEMinusInvPin < 0.05;
1553 +      }
1554 +  }
1555 +
1556 +  else if(variable == "genMatchedId"){
1557 +    int index = getGenMatchedParticleIndex(object);
1558 +    if(index == -1) value = 0;
1559 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1560 +  }
1561 +  else if(variable == "genMatchedMotherId"){
1562 +    int index = getGenMatchedParticleIndex(object);
1563 +    if(index == -1) value = 0;
1564 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1565 +  }
1566 +  else if(variable == "genMatchedGrandmotherId"){
1567 +    int index = getGenMatchedParticleIndex(object);
1568 +    if(index == -1) value = 0;
1569 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1570 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1571 +      if(motherIndex == -1) value = 0;
1572 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1573 +    }
1574 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1575 +  }
1576 +
1577 +
1578  
1579    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1580 <  
1580 >
1581 >  value = applyFunction(function, value);
1582  
1583    return value;
1584   }
# Line 933 | Line 1653 | OSUAnalysis::valueLookup (const BNevent*
1653    else if(variable == "id1") value = object->id1;
1654    else if(variable == "id2") value = object->id2;
1655    else if(variable == "evt") value = object->evt;
1656 +  else if(variable == "puScaleFactor"){
1657 +    if(datasetType_ != "data")
1658 +      value = puWeight_->at (events->at (0).numTruePV);
1659 +    else
1660 +      value = 1.0;
1661 +  }
1662 +  else if(variable == "muonScaleFactor"){
1663 +    if(datasetType_ != "data")
1664 +      value = muonSFWeight_->at (chosenMuon ()->eta);
1665 +    else
1666 +      value = 1.0;
1667 +  }
1668 +  else if(variable == "electronScaleFactor"){
1669 +    if(datasetType_ != "data")
1670 +      value = electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
1671 +    else
1672 +      value = 1.0;
1673 +  }
1674  
1675    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1676 <  
1676 >
1677 >  value = applyFunction(function, value);
1678  
1679    return value;
1680   }
# Line 986 | Line 1725 | OSUAnalysis::valueLookup (const BNtau* o
1725    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1726  
1727  
1728 +  else if(variable == "genMatchedId"){
1729 +    int index = getGenMatchedParticleIndex(object);
1730 +    if(index == -1) value = 0;
1731 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1732 +  }
1733 +  else if(variable == "genMatchedMotherId"){
1734 +    int index = getGenMatchedParticleIndex(object);
1735 +    if(index == -1) value = 0;
1736 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1737 +  }
1738 +  else if(variable == "genMatchedGrandmotherId"){
1739 +    int index = getGenMatchedParticleIndex(object);
1740 +    if(index == -1) value = 0;
1741 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1742 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1743 +      if(motherIndex == -1) value = 0;
1744 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1745 +    }
1746 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1747 +  }
1748 +
1749 +
1750    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1751 <  
1751 >
1752 >  value = applyFunction(function, value);
1753  
1754    return value;
1755   }
# Line 1058 | Line 1820 | OSUAnalysis::valueLookup (const BNmet* o
1820    else if(variable == "pfT1jet10phi") value = object->pfT1jet10phi;
1821  
1822    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1823 <  
1823 >
1824 >  value = applyFunction(function, value);
1825  
1826    return value;
1827   }
# Line 1067 | Line 1830 | double
1830   OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){
1831  
1832    double value = 0.0;
1833 +  double pMag = sqrt(object->pt * object->pt +
1834 +                         object->pz * object->pz);
1835  
1836    if(variable == "pt") value = object->pt;
1837    else if(variable == "px") value = object->px;
# Line 1087 | Line 1852 | OSUAnalysis::valueLookup (const BNtrack*
1852    else if(variable == "isHighPurity") value = object->isHighPurity;
1853  
1854  
1855 <  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1855 >  //additional BNs info for disappTrks
1856 >  else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
1857 >  else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
1858 >  else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
1859 >  else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
1860 >  else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
1861 >  else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
1862 >  else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
1863 >  else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
1864 >  else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
1865 >  else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
1866    
1867  
1868 +  //user defined variables
1869 +  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;
1870 +  else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
1871 +  else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
1872 +  else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
1873 +
1874 +
1875 +
1876 +
1877 +
1878 +  else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
1879 +  else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
1880 +  else if(variable == "ptErrorByPt") value = (object->ptError/object->pt);
1881 +  else if(variable == "ptError") value = object->ptError;
1882 +  else if(variable == "ptRes") value = getTrkPtRes(object);
1883 +
1884 +  else if(variable == "genDeltaRLowest") value = getGenDeltaRLowest(object);
1885 +
1886 +  else if(variable == "genId"){
1887 +    int index = getGenMatchedParticleIndex(object);
1888 +    if(index == -1) value = 0;
1889 +    else value = mcparticles->at(index).id;
1890 +  }
1891 +
1892 +
1893 +  else if(variable == "genMatchedId"){
1894 +    int index = getGenMatchedParticleIndex(object);
1895 +    if(index == -1) value = 0;
1896 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1897 +  }
1898 +  else if(variable == "genMatchedMotherId"){
1899 +    int index = getGenMatchedParticleIndex(object);
1900 +    if(index == -1) value = 0;
1901 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1902 +  }
1903 +  else if(variable == "genMatchedGrandmotherId"){
1904 +    int index = getGenMatchedParticleIndex(object);
1905 +    if(index == -1) value = 0;
1906 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1907 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1908 +      if(motherIndex == -1) value = 0;
1909 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1910 +    }
1911 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1912 +  }
1913 +
1914 +
1915 +
1916 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1917 +
1918 +  value = applyFunction(function, value);
1919 +
1920    return value;
1921   }
1922  
# Line 1115 | Line 1942 | OSUAnalysis::valueLookup (const BNgenjet
1942  
1943  
1944    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1945 <  
1945 >
1946 >  value = applyFunction(function, value);
1947  
1948    return value;
1949   }
# Line 1210 | Line 2038 | OSUAnalysis::valueLookup (const BNmcpart
2038    else if(variable == "grandMother11Status") value = object->grandMother11Status;
2039    else if(variable == "grandMother11Charge") value = object->grandMother11Charge;
2040  
2041 +  //user-defined variables
2042 +  else if (variable == "d0"){
2043 +    double vx = object->vx - chosenVertex ()->x,
2044 +      vy = object->vy - chosenVertex ()->y,
2045 +      px = object->px,
2046 +      py = object->py,
2047 +      pt = object->pt;
2048 +    value = (-vx * py + vy * px) / pt;
2049 +  }
2050 +  else if (variable == "dz"){
2051 +    double vx = object->vx - chosenVertex ()->x,
2052 +      vy = object->vy - chosenVertex ()->y,
2053 +      vz = object->vz - chosenVertex ()->z,
2054 +      px = object->px,
2055 +      py = object->py,
2056 +      pz = object->pz,
2057 +      pt = object->pt;
2058 +    value = vz - (vx * px + vy * py)/pt * (pz/pt);
2059 +  }
2060 +  else if(variable == "v0"){
2061 +    value = sqrt(object->vx*object->vx + object->vy*object->vy);
2062 +  }
2063 +  else if(variable == "deltaV0"){
2064 +    value = sqrt((object->vx-chosenVertex ()->x)*(object->vx-chosenVertex ()->x) + (object->vy-chosenVertex ()->y)*(object->vy-chosenVertex ()->y));
2065 +  }
2066 +  else if (variable == "deltaVx"){
2067 +    value = object->vx - chosenVertex ()->x;
2068 +  }
2069 +  else if (variable == "deltaVy"){
2070 +    value = object->vy - chosenVertex ()->y;
2071 +  }
2072 +  else if (variable == "deltaVz"){
2073 +    value = object->vz - chosenVertex ()->z;
2074 +  }
2075 +
2076 +
2077    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2078 <  
2078 >
2079 >  value = applyFunction(function, value);
2080  
2081    return value;
2082   }
# Line 1237 | Line 2102 | OSUAnalysis::valueLookup (const BNprimar
2102  
2103  
2104    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2105 <  
2105 >
2106 >  value = applyFunction(function, value);
2107  
2108    return value;
2109   }
# Line 1253 | Line 2119 | OSUAnalysis::valueLookup (const BNbxlumi
2119  
2120  
2121    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2122 <  
2122 >
2123 >  value = applyFunction(function, value);
2124  
2125    return value;
2126   }
# Line 1334 | Line 2201 | OSUAnalysis::valueLookup (const BNphoton
2201    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
2202  
2203  
2204 +
2205 +
2206 +
2207 +
2208 +
2209 +  else if(variable == "genMatchedId"){
2210 +    int index = getGenMatchedParticleIndex(object);
2211 +    if(index == -1) value = 0;
2212 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
2213 +  }
2214 +  else if(variable == "genMatchedMotherId"){
2215 +    int index = getGenMatchedParticleIndex(object);
2216 +    if(index == -1) value = 0;
2217 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2218 +  }
2219 +  else if(variable == "genMatchedGrandmotherId"){
2220 +    int index = getGenMatchedParticleIndex(object);
2221 +    if(index == -1) value = 0;
2222 +    else if(fabs(mcparticles->at(index).motherId) == 15){
2223 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2224 +      if(motherIndex == -1) value = 0;
2225 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2226 +    }
2227 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2228 +  }
2229 +
2230 +
2231    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2232 <  
2232 >
2233 >  value = applyFunction(function, value);
2234  
2235    return value;
2236   }
# Line 1356 | Line 2251 | OSUAnalysis::valueLookup (const BNsuperc
2251  
2252  
2253    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2254 <  
2254 >
2255 >  value = applyFunction(function, value);
2256 >
2257 >  return value;
2258 > }
2259 >
2260 >
2261 > double
2262 > OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function){
2263 >
2264 >  double value = 0.0;
2265 >
2266 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2267 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2268 >  else if(variable == "invMass"){
2269 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2270 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2271 >    value = (fourVector1 + fourVector2).M();
2272 >  }
2273 >  else if(variable == "threeDAngle")
2274 >    {
2275 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2276 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2277 >      value = (threeVector1.Angle(threeVector2));
2278 >    }
2279 >  else if(variable == "alpha")
2280 >    {
2281 >      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
2282 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2283 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2284 >      value = (pi-threeVector1.Angle(threeVector2));
2285 >    }
2286 >  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2287 >  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2288 >  else if(variable == "d0Sign"){
2289 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2290 >    if(d0Sign < 0) value = -0.5;
2291 >    else if (d0Sign > 0) value = 0.5;
2292 >    else value = -999;
2293 >  }
2294 >  else if(variable == "chargeProduct"){
2295 >    value = object1->charge*object2->charge;
2296 >  }
2297 >  else if(variable == "muon1CorrectedD0Vertex"){
2298 >    value = object1->correctedD0Vertex;
2299 >  }
2300 >  else if(variable == "muon2CorrectedD0Vertex"){
2301 >    value = object2->correctedD0Vertex;
2302 >  }
2303 > else if(variable == "muon1timeAtIpInOut"){
2304 >    value = object1->timeAtIpInOut;
2305 >  }
2306 > else if(variable == "muon2timeAtIpInOut"){
2307 >    value = object2->timeAtIpInOut;
2308 >  }
2309 > else if(variable == "muon1correctedD0")
2310 >   {
2311 >     value = object1->correctedD0;
2312 >   }
2313 > else if(variable == "muon1correctedD0Vertex")
2314 >   {
2315 >     value = object1->correctedD0Vertex;
2316 >   }
2317 >  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2318 >
2319 >  value = applyFunction(function, value);
2320  
2321    return value;
2322   }
2323  
2324 + double
2325 + OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function){
2326 +
2327 +  double value = 0.0;
2328 +
2329 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2330 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2331 +  else if(variable == "invMass"){
2332 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2333 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2334 +    value = (fourVector1 + fourVector2).M();
2335 +  }
2336 +  else if(variable == "threeDAngle")
2337 +    {
2338 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2339 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2340 +      value = (threeVector1.Angle(threeVector2));
2341 +    }
2342 +  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2343 +  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2344 +  else if(variable == "d0Sign"){
2345 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2346 +    if(d0Sign < 0) value = -0.5;
2347 +    else if (d0Sign > 0) value = 0.5;
2348 +    else value = -999;
2349 +  }
2350 +  else if(variable == "chargeProduct"){
2351 +    value = object1->charge*object2->charge;
2352 +  }
2353 +  else if(variable == "electron1CorrectedD0Vertex"){
2354 +    value = object1->correctedD0Vertex;
2355 +  }
2356 +  else if(variable == "electron2CorrectedD0Vertex"){
2357 +    value = object2->correctedD0Vertex;
2358 +  }
2359 +
2360 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2361 +
2362 +  value = applyFunction(function, value);
2363  
2364 < template <class InputCollection>
2364 >  return value;
2365 > }
2366 >
2367 > double
2368 > OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function){
2369 >
2370 >  double value = 0.0;
2371 >
2372 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2373 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2374 >  else if(variable == "invMass"){
2375 >    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2376 >    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2377 >    value = (fourVector1 + fourVector2).M();
2378 >  }
2379 >  else if(variable == "threeDAngle")
2380 >    {
2381 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2382 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2383 >      value = (threeVector1.Angle(threeVector2));
2384 >    }
2385 >  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2386 >  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2387 >  else if(variable == "d0Sign"){
2388 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2389 >    if(d0Sign < 0) value = -0.5;
2390 >    else if (d0Sign > 0) value = 0.5;
2391 >    else value = -999;
2392 >  }
2393 >  else if(variable == "chargeProduct"){
2394 >    value = object1->charge*object2->charge;
2395 >  }
2396 >  else if(variable == "electronCorrectedD0Vertex"){
2397 >    value = object1->correctedD0Vertex;
2398 >  }
2399 >  else if(variable == "muonCorrectedD0Vertex"){
2400 >    value = object2->correctedD0Vertex;
2401 >  }
2402 >  else if(variable == "electronCorrectedD0"){
2403 >    value = object1->correctedD0;
2404 >  }
2405 >  else if(variable == "muonCorrectedD0"){
2406 >    value = object2->correctedD0;
2407 >  }
2408 >  else if(variable == "electronDetIso"){
2409 >    value = (object1->trackIso) / object1->pt;
2410 >  }
2411 >  else if(variable == "muonDetIso"){
2412 >    value = (object2->trackIsoDR03) / object2->pt;
2413 >  }
2414 >  else if(variable == "chargeProduct"){
2415 >    value = object1->charge * object2->charge;
2416 >  }
2417 >
2418 >
2419 >  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2420 >
2421 >  value = applyFunction(function, value);
2422 >
2423 >  return value;
2424 > }
2425 >
2426 >
2427 > // Calculate the number of tracks in cone of DeltaR<0.5 around track1.
2428 > // Return true iff no other tracks are found in this cone.
2429 > int
2430 > OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
2431 >  for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
2432 >    if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.
2433 >    double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
2434 >    if(deltaRtrk < 0.5) return 0;
2435 >  }
2436 >  return 1;
2437 >
2438 > }
2439 >
2440 >
2441 > double
2442 > OSUAnalysis::getTrkPtRes (const BNtrack* track1){
2443 >
2444 >  double ptTrue = getTrkPtTrue(track1, mcparticles.product());
2445 >  double PtRes = (track1->pt - ptTrue) / ptTrue;
2446 >
2447 >  return PtRes;
2448 >
2449 > }
2450 >
2451 >
2452 > double
2453 > OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
2454 >  double value = -99;
2455 >  double genDeltaRLowest = 999;
2456 >
2457 >  for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
2458 >    double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
2459 >    if (genDeltaRtemp < genDeltaRLowest) {
2460 >      genDeltaRLowest = genDeltaRtemp;
2461 >      if (genDeltaRLowest < 0.05) {   // Only consider it truth-matched if DeltaR<0.15.
2462 >        double ptTrue = genPart->pt;
2463 >        value = ptTrue;
2464 >      }
2465 >    }
2466 >  }
2467 >
2468 >  return value;
2469 >
2470 > }
2471 >
2472 > //creates a map of the dead Ecal channels in the barrel and endcap
2473 > //to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
2474 > //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
2475 > void
2476 > OSUAnalysis::WriteDeadEcal (){
2477 >  double etaEcal, phiEcal;
2478 >  ifstream DeadEcalFile(deadEcalFile_);
2479 >  if(!DeadEcalFile) {
2480 >    cout << "Error: DeadEcalFile has not been found." << endl;
2481 >    return;
2482 >  }
2483 >  if(DeadEcalVec.size()!= 0){
2484 >    cout << "Error: DeadEcalVec has a nonzero size" << endl;
2485 >    return;
2486 >  }
2487 >  while(!DeadEcalFile.eof())
2488 >    {
2489 >      DeadEcalFile >> etaEcal >> phiEcal;
2490 >      DeadEcal newChan;
2491 >      newChan.etaEcal = etaEcal;
2492 >      newChan.phiEcal = phiEcal;
2493 >      DeadEcalVec.push_back(newChan);
2494 >    }
2495 >  if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl;
2496 > }
2497 >
2498 > //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
2499 > int
2500 > OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
2501 >  double deltaRLowest = 999;
2502 >  int value = 0;
2503 >  if (DeadEcalVec.size() == 0) WriteDeadEcal();
2504 >  for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
2505 >    double eta = ecal->etaEcal;
2506 >    double phi = ecal->phiEcal;
2507 >    double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
2508 >    if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
2509 >  }
2510 >  if (deltaRLowest<0.05) {value = 1;}
2511 >  else {value = 0;}
2512 >  return value;
2513 > }
2514 >
2515 > // Returns the smallest DeltaR between the object and any generated true particle in the event.  
2516 > template <class InputObject>
2517 > double OSUAnalysis::getGenDeltaRLowest(InputObject object){
2518 >  double genDeltaRLowest = 999.;
2519 >  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2520 >    double deltaRtemp = deltaR(mcparticle->eta, mcparticle->phi, object->eta, object->phi);
2521 >    if (deltaRtemp < genDeltaRLowest) genDeltaRLowest = deltaRtemp;
2522 > }
2523 >  return genDeltaRLowest;
2524 > }
2525 >
2526 > double
2527 > OSUAnalysis::applyFunction(string function, double value){
2528 >
2529 >  if(function == "abs") value = fabs(value);
2530 >  else if(function == "fabs") value = fabs(value);
2531 >  else if(function == "log10") value = log10(value);
2532 >  else if(function == "log") value = log10(value);
2533 >
2534 >  else if(function == "") value = value;
2535 >  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
2536 >
2537 >  return value;
2538 >
2539 > }
2540 >
2541 >
2542 > template <class InputCollection>
2543   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2544  
2545  
2546    for (uint object = 0; object != inputCollection->size(); object++){
2547 <    
2547 >
2548 >
2549      bool decision = true;//object passes if this cut doesn't cut on that type of object
1372    
1373    if(currentCut.inputCollection == inputType){
1374      
1375      double value = valueLookup(&inputCollection->at(object), currentCut.variable, currentCut.function);
2550  
2551 <      decision = evaluateComparison(value,currentCut.comparativeOperator,currentCut.cutValue);
2551 >
2552 >    if(currentCut.inputCollection == inputType){
2553 >
2554 >      vector<bool> subcutDecisions;
2555 >      for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2556 >        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2557 >        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2558 >      }
2559 >      if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2560 >      else{
2561 >        bool tempDecision = true;
2562 >        for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2563 >          if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2564 >            tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2565 >          else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2566 >            tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
2567 >        }
2568 >        decision = tempDecision;
2569 >      }
2570      }
2571 <    individualFlags[inputType].at(currentCutIndex).push_back(decision);
2572 <    
2571 >    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2572 >
2573  
2574      //set flags for objects that pass each cut AND all the previous cuts
2575      bool previousCumulativeFlag = true;
2576      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2577 <      if(previousCumulativeFlag && individualFlags[inputType].at(previousCutIndex).at(object)) previousCumulativeFlag = true;
2577 >      if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(object)) previousCumulativeFlag = true;
2578        else{ previousCumulativeFlag = false; break;}
2579      }
2580 <    cumulativeFlags[inputType].at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2581 <    
2580 >    cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2581 >
2582 >
2583    }
2584 <  
2584 >
2585   }
2586  
2587  
2588 + template <class InputCollection1, class InputCollection2>
2589 + void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
2590 +                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
2591  
2592  
2593 +  bool sameObjects = false;
2594 +  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2595 +
2596 +
2597 +  int counter = 0;
2598 +  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2599 +    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2600 +
2601 +      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2602  
1398 DEFINE_FWK_MODULE(OSUAnalysis);
2603  
2604 +      bool decision = true;//object passes if this cut doesn't cut on that type of object
2605  
2606 +      if(currentCut.inputCollection == inputType){
2607 +
2608 +        vector<bool> subcutDecisions;
2609 +        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2610 +          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2611 +          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2612 +        }
2613 +
2614 +        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2615 +        else{
2616 +          bool tempDecision = subcutDecisions.at(0);
2617 +          for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
2618 +            if(currentCut.logicalOperators.at(subcutIndex-1) == "&" || currentCut.logicalOperators.at(subcutIndex-1) == "&&")
2619 +              tempDecision = tempDecision && subcutDecisions.at(subcutIndex);
2620 +            else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
2621 +              tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
2622 +          }
2623 +          decision = tempDecision;
2624 +        }
2625 +      }
2626 +      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2627 +
2628 +      //set flags for objects that pass each cut AND all the previous cuts
2629 +      bool previousCumulativeFlag = true;
2630 +      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2631 +        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
2632 +        else{ previousCumulativeFlag = false; break;}
2633 +      }
2634 +      //apply flags for the components of the composite object as well
2635 +      bool currentCumulativeFlag = true;
2636 +      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
2637 +      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
2638 +      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
2639 +      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
2640 +      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
2641 +
2642 +      counter++;
2643 +    }
2644 +
2645 +  }
2646 +
2647 +
2648 + }
2649 +
2650 +
2651 + template <class InputCollection>
2652 + void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
2653 +
2654 +  for (uint object = 0; object != inputCollection->size(); object++){
2655 +
2656 +    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2657 +
2658 +    string currentString = parameters.inputVariables.at(0);
2659 +    string inputVariable = "";
2660 +    string function = "";
2661 +    if(currentString.find("(")==std::string::npos){
2662 +      inputVariable = currentString;// variable to cut on
2663 +    }
2664 +    else{
2665 +      function = currentString.substr(0,currentString.find("("));//function comes before the "("
2666 +      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2667 +      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2668 +    }
2669 +
2670 +    double value = valueLookup(&inputCollection->at(object), inputVariable, function);
2671 +    histo->Fill(value,scaleFactor);
2672 +
2673 +  }
2674 +
2675 + }
2676 +
2677 + template <class InputCollection1, class InputCollection2>
2678 + void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
2679 +
2680 +  bool sameObjects = false;
2681 +  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2682 +
2683 +  int pairCounter = 0;
2684 +  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2685 +    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2686 +
2687 +      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2688 +
2689 +      //only take objects which have passed all cuts and pairs which have passed all cuts
2690 +      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2691 +      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2692 +      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2693 +
2694 +      string currentString = parameters.inputVariables.at(0);
2695 +      string inputVariable = "";
2696 +      string function = "";
2697 +      if(currentString.find("(")==std::string::npos){
2698 +        inputVariable = currentString;// variable to cut on
2699 +      }
2700 +      else{
2701 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("
2702 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2703 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2704 +      }
2705 +
2706 +      double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2707 +      histo->Fill(value,scaleFactor);
2708 +
2709 +      pairCounter++;
2710 +    }
2711 +  }
2712 +
2713 + }
2714 +
2715 +
2716 + template <class InputCollection>
2717 + void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
2718 +
2719 +  for (uint object = 0; object != inputCollection->size(); object++){
2720 +
2721 +    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2722 +
2723 +    string currentString = parameters.inputVariables.at(0);
2724 +    string inputVariable = "";
2725 +    string function = "";
2726 +    if(currentString.find("(")==std::string::npos){
2727 +      inputVariable = currentString;// variable to cut on
2728 +    }
2729 +    else{
2730 +      function = currentString.substr(0,currentString.find("("));//function comes before the "("
2731 +      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2732 +      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2733 +    }
2734 +    double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
2735 +
2736 +    currentString = parameters.inputVariables.at(1);
2737 +    inputVariable = "";
2738 +    function = "";
2739 +    if(currentString.find("(")==std::string::npos){
2740 +      inputVariable = currentString;// variable to cut on
2741 +    }
2742 +    else{
2743 +      function = currentString.substr(0,currentString.find("("));//function comes before the "("
2744 +      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2745 +      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2746 +    }
2747 +
2748 +    double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
2749 +
2750 +    histo->Fill(valueX,valueY,scaleFactor);
2751 +
2752 +  }
2753 +
2754 + }
2755 +
2756 + template <class InputCollection1, class InputCollection2>
2757 + void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
2758 +
2759 +  bool sameObjects = false;
2760 +  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2761 +
2762 +  int pairCounter = 0;
2763 +  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2764 +    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2765 +
2766 +      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2767 +
2768 +      //only take objects which have passed all cuts and pairs which have passed all cuts
2769 +      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2770 +      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2771 +      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2772 +
2773 +      string currentString = parameters.inputVariables.at(0);
2774 +      string inputVariable = "";
2775 +      string function = "";
2776 +      if(currentString.find("(")==std::string::npos){
2777 +        inputVariable = currentString;// variable to cut on
2778 +      }
2779 +      else{
2780 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("
2781 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2782 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2783 +      }
2784 +      double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2785 +
2786 +      currentString = parameters.inputVariables.at(1);
2787 +      inputVariable = "";
2788 +      function = "";
2789 +      if(currentString.find("(")==std::string::npos){
2790 +        inputVariable = currentString;// variable to cut on
2791 +      }
2792 +      else{
2793 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("
2794 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
2795 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
2796 +      }
2797 +      double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2798 +
2799 +
2800 +      histo->Fill(valueX,valueY,scaleFactor);
2801 +
2802 +      pairCounter++;
2803 +
2804 +    }
2805 +  }
2806 +
2807 + }
2808 +
2809 +
2810 + template <class InputObject>
2811 + int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
2812 +
2813 +  int bestMatchIndex = -1;
2814 +  double bestMatchDeltaR = 999;
2815 +
2816 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2817 +
2818 +    double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
2819 +    if(currentDeltaR > 0.05) continue;
2820 + //     cout << std::setprecision(3) << std::setw(20)
2821 + //          << "\tcurrentParticle:  eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
2822 + //          << std::setw(20)
2823 + //          << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
2824 + //          << std::setw(20)
2825 + //          << "\tdeltaR = " << currentDeltaR
2826 + //          << std::setprecision(1)
2827 + //          << std::setw(20)
2828 + //          << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
2829 + //          << std::setw(20)
2830 + //          << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
2831 + //          << std::setw(20)
2832 + //          << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
2833 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2834 +      bestMatchIndex = mcparticle - mcparticles->begin();
2835 +      bestMatchDeltaR = currentDeltaR;
2836 +    }
2837 +
2838 +  }
2839 + //   if(bestMatchDeltaR != 999)  cout << "bestMatch:  deltaR = " << bestMatchDeltaR << "   id = " << mcparticles->at(bestMatchIndex).id << "   motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
2840 + //   else cout << "no match found..." << endl;
2841 +  return bestMatchIndex;
2842 +
2843 + }
2844 +
2845 +
2846 + int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
2847 +
2848 +  int bestMatchIndex = -1;
2849 +  double bestMatchDeltaR = 999;
2850 +
2851 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2852 +
2853 +    if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
2854 +
2855 +    double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
2856 +    if(currentDeltaR > 0.05) continue;
2857 +
2858 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2859 +      bestMatchIndex = mcparticle - mcparticles->begin();
2860 +      bestMatchDeltaR = currentDeltaR;
2861 +    }
2862 +
2863 +  }
2864 +
2865 +  return bestMatchIndex;
2866 + }
2867 +
2868 +
2869 + // bin      particle type
2870 + // ---      -------------
2871 + //  0        unmatched
2872 + //  1        u
2873 + //  2        d
2874 + //  3        s
2875 + //  4        c
2876 + //  5        b
2877 + //  6        t
2878 + //  7        e
2879 + //  8        mu
2880 + //  9        tau
2881 + // 10        nu
2882 + // 11        g
2883 + // 12        gamma
2884 + // 13        Z
2885 + // 14        W
2886 + // 15        light meson
2887 + // 16        K meson
2888 + // 17        D meson
2889 + // 18        B meson
2890 + // 19        light baryon
2891 + // 20        strange baryon
2892 + // 21        charm baryon
2893 + // 22        bottom baryon
2894 + // 23        other
2895 +
2896 + int OSUAnalysis::getPdgIdBinValue(int pdgId){
2897 +
2898 +  int binValue = -999;
2899 +  int absPdgId = fabs(pdgId);
2900 +  if(pdgId == -1) binValue = 0;
2901 +  else if(absPdgId <= 6 ) binValue = absPdgId;
2902 +  else if(absPdgId == 11 ) binValue = 7;
2903 +  else if(absPdgId == 13 ) binValue = 8;
2904 +  else if(absPdgId == 15 ) binValue = 9;
2905 +  else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
2906 +  else if(absPdgId == 21 ) binValue = 11;
2907 +  else if(absPdgId == 22 ) binValue = 12;
2908 +  else if(absPdgId == 23 ) binValue = 13;
2909 +  else if(absPdgId == 24 ) binValue = 14;
2910 +  else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130  ) binValue = 15;
2911 +  else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400)  ) binValue = 16;
2912 +  else if(absPdgId > 400 && absPdgId < 500  ) binValue = 17;
2913 +  else if(absPdgId > 500 && absPdgId < 600  ) binValue = 18;
2914 +  else if(absPdgId > 1000 && absPdgId < 3000  ) binValue = 19;
2915 +  else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
2916 +  else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
2917 +  else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
2918 +
2919 +  else binValue = 23;
2920 +
2921 +  return binValue;
2922 +
2923 + }
2924 +
2925 + const BNprimaryvertex *
2926 + OSUAnalysis::chosenVertex ()
2927 + {
2928 +  const BNprimaryvertex *chosenVertex = 0;
2929 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
2930 +    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back().size() ? cumulativeFlags.at("primaryvertexs").back() :
2931 +                               cumulativeFlags.at("primaryvertexs").at(cumulativeFlags.at("primaryvertexs").size() - 2);
2932 +    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
2933 +      if(!vertexFlags.at(vertexIndex)) continue;
2934 +      chosenVertex = & primaryvertexs->at(vertexIndex);
2935 +      break;
2936 +    }
2937 +  }
2938 +  else {
2939 +    chosenVertex = &primaryvertexs->at(0);
2940 +  }
2941 +
2942 +  return chosenVertex;
2943 + }
2944 +
2945 + const BNmet *
2946 + OSUAnalysis::chosenMET ()
2947 + {
2948 +  const BNmet *chosenMET = 0;
2949 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "mets") != objectsToCut.end()) {
2950 +    vector<bool> metFlags = cumulativeFlags.at("mets").back().size() ? cumulativeFlags.at("mets").back() :
2951 +                            cumulativeFlags.at("mets").at(cumulativeFlags.at("mets").size() - 2);
2952 +    for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
2953 +      if(!metFlags.at(metIndex)) continue;
2954 +      chosenMET = & mets->at(metIndex);
2955 +      break;
2956 +    }
2957 +  }
2958 +  else {
2959 +    chosenMET = &mets->at(0);
2960 +  }
2961 +
2962 +  return chosenMET;
2963 + }
2964 +
2965 + const BNelectron *
2966 + OSUAnalysis::chosenElectron ()
2967 + {
2968 +  const BNelectron *chosenElectron = 0;
2969 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "electrons") != objectsToCut.end()) {
2970 +    vector<bool> electronFlags = cumulativeFlags.at("electrons").back().size() ? cumulativeFlags.at("electrons").back() :
2971 +                                 cumulativeFlags.at("electrons").at(cumulativeFlags.at("electrons").size() - 2);
2972 +    for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
2973 +      if(!electronFlags.at(electronIndex)) continue;
2974 +      chosenElectron = & electrons->at(electronIndex);
2975 +      break;
2976 +    }
2977 +  }
2978 +  else {
2979 +    chosenElectron = &electrons->at(0);
2980 +  }
2981 +
2982 +  return chosenElectron;
2983 + }
2984 +
2985 + const BNmuon *
2986 + OSUAnalysis::chosenMuon ()
2987 + {
2988 +  const BNmuon *chosenMuon = 0;
2989 +  if(std::find(objectsToCut.begin(), objectsToCut.end(), "muons") != objectsToCut.end()) {
2990 +    vector<bool> muonFlags = cumulativeFlags.at("muons").back().size() ? cumulativeFlags.at("muons").back() :
2991 +                             cumulativeFlags.at("muons").at(cumulativeFlags.at("muons").size() - 2);
2992 +    for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
2993 +      if(!muonFlags.at(muonIndex)) continue;
2994 +      chosenMuon = & muons->at(muonIndex);
2995 +      break;
2996 +    }
2997 +  }
2998 +  else {
2999 +    chosenMuon = &muons->at(0);
3000 +  }
3001 +
3002 +  return chosenMuon;
3003 + }
3004 +
3005 +
3006 + DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines