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.1 by lantonel, Mon Jan 21 10:39:34 2013 UTC vs.
Revision 1.28 by jbrinson, Thu Mar 14 10:06:38 2013 UTC

# Line 15 | Line 15 | OSUAnalysis::OSUAnalysis (const edm::Par
15    bxlumis_ (cfg.getParameter<edm::InputTag> ("bxlumis")),
16    photons_ (cfg.getParameter<edm::InputTag> ("photons")),
17    superclusters_ (cfg.getParameter<edm::InputTag> ("superclusters")),
18 <
19 <
20 <  channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels"))
18 >  triggers_ (cfg.getParameter<edm::InputTag> ("triggers")),
19 >  puFile_ (cfg.getParameter<std::string> ("puFile")),
20 >  deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")),
21 >  dataPU_ (cfg.getParameter<std::string> ("dataPU")),
22 >  dataset_ (cfg.getParameter<std::string> ("dataset")),
23 >  datasetType_ (cfg.getParameter<std::string> ("datasetType")),
24 >  channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
25 >  histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
26 >  plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents"))
27  
28   {
29 +
30    TH1::SetDefaultSumw2 ();
31  
32 +  //create pile-up reweighting object, if necessary
33 +  //  if(datasetType_ != "data") puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
34 +
35    // Construct Cutflow Objects. These store the results of cut decisions and
36    // handle filling cut flow histograms.
37    masterCutFlow_ = new CutFlow (fs_);
38    std::vector<TFileDirectory> directories;
39  
40 +  //always get vertex collection so we can assign the primary vertex in the event
41 +  
42 + objectsToGet.push_back("primaryvertexs");
43 +  //always make the plot of number of primary verticex (to check pile-up reweighting)
44 +  objectsToPlot.push_back("primaryvertexs");
45 +
46 +  //always get the MC particles to do GEN-matching
47 +  objectsToGet.push_back("mcparticles");
48 +
49 +  //always get the event collection to do pile-up reweighting
50 +  objectsToGet.push_back("events");
51 +
52 +  //parse the histogram definitions
53 +  for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
54 +
55 +    string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
56 +    if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
57 +    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
58 +      objectsToGet.push_back(tempInputCollection);
59 +      objectsToPlot.push_back(tempInputCollection);
60 +      objectsToCut.push_back(tempInputCollection);
61 +    }
62 +    else{//pair of objects, need to add them both to the things to objectsToGet
63 +      int dashIndex = tempInputCollection.find("-");
64 +      int spaceIndex = tempInputCollection.find(" ");
65 +      int secondWordLength = spaceIndex - dashIndex;
66 +      objectsToGet.push_back(tempInputCollection);
67 +      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
68 +      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
69 +      objectsToPlot.push_back(tempInputCollection);
70 +      objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
71 +      objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
72 +      objectsToCut.push_back(tempInputCollection);
73 +      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
74 +      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
75 +
76 +      }
77 +
78 +    vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
79 +    
80 +    for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
81 +
82 +      histogram tempHistogram;
83 +      tempHistogram.inputCollection = tempInputCollection;
84 +      tempHistogram.name = histogramList_.at(currentHistogram).getParameter<string>("name");
85 +      tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
86 +      tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
87 +      tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
88 +      
89 +      histograms.push_back(tempHistogram);
90 +
91 +    }
92 +  }
93 +  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
94 +  sort( objectsToPlot.begin(), objectsToPlot.end() );
95 +  objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
96 +
97 +
98 +
99 +  //add histograms with the gen-matched id, mother id, and grandmother id
100 +  for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
101 +
102 +    string currentObject = objectsToPlot.at(currentObjectIndex);
103 +    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
104 +
105 +    histogram tempIdHisto;
106 +    histogram tempMomIdHisto;
107 +    histogram tempGmaIdHisto;
108 +
109 +    tempIdHisto.inputCollection = currentObject;
110 +    tempMomIdHisto.inputCollection = currentObject;
111 +    tempGmaIdHisto.inputCollection = currentObject;
112 +
113 +    currentObject = currentObject.substr(0, currentObject.size()-1);
114 +    tempIdHisto.name = currentObject+"GenMatchId";
115 +    tempMomIdHisto.name = currentObject+"GenMatchMotherId";
116 +    tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
117 +
118 +    currentObject.at(0) = toupper(currentObject.at(0));
119 +    tempIdHisto.title = currentObject+" Gen-matched Particle";
120 +    tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
121 +    tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
122 +
123 +    int maxNum = 24;
124 +    vector<double> binVector;
125 +    binVector.push_back(maxNum);
126 +    binVector.push_back(0);
127 +    binVector.push_back(maxNum);
128 +
129 +    tempIdHisto.bins = binVector;
130 +    tempIdHisto.inputVariables.push_back("genMatchedId");
131 +    tempMomIdHisto.bins = binVector;
132 +    tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
133 +    tempGmaIdHisto.bins = binVector;
134 +    tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
135 +
136 +    histograms.push_back(tempIdHisto);
137 +    histograms.push_back(tempMomIdHisto);
138 +    histograms.push_back(tempGmaIdHisto);
139 +    
140 +  }
141 +
142 +
143  
144    channel tempChannel;
145    //loop over all channels (event selections)
146    for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
147 <    
147 >
148      //get name of channel
149      string channelName  (channels_.at(currentChannel).getParameter<string>("name"));
150      tempChannel.name = channelName;
151      TString channelLabel = channelName;
152  
153 +    //set triggers for this channel
154 +    vector<string> triggerNames;
155 +    triggerNames.clear();
156 +    tempChannel.triggers.clear();
157 +    if(channels_.at(currentChannel).exists("triggers")){
158 +      triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
159 +      for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
160 +        tempChannel.triggers.push_back(triggerNames.at(trigger));
161 +      objectsToGet.push_back("triggers");
162 +    }
163 +
164 +
165 +
166  
167      //create cutFlow for this channel
168      cutFlows_.push_back (new CutFlow (fs_, channelName));
# Line 44 | Line 170 | OSUAnalysis::OSUAnalysis (const edm::Par
170      //book a directory in the output file with the name of the channel
171      TFileDirectory subDir = fs_->mkdir( channelName );
172      directories.push_back(subDir);
47    std::map<std::string, TH1D*> histoMap;
48    oneDHists_.push_back(histoMap);
173  
174 +    std::map<std::string, TH1D*> oneDhistoMap;
175 +    oneDHists_.push_back(oneDhistoMap);
176 +    std::map<std::string, TH2D*> twoDhistoMap;
177 +    twoDHists_.push_back(twoDhistoMap);
178  
51    //muon histograms
52    oneDHists_.at(currentChannel)["muonPt"]  = directories.at(currentChannel).make<TH1D> ("muonPt",channelLabel+" channel: Muon Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
53    oneDHists_.at(currentChannel)["muonEta"] = directories.at(currentChannel).make<TH1D> ("muonEta",channelLabel+" channel: Muon Eta; #eta", 100, -5, 5);
54    oneDHists_.at(currentChannel)["muonD0"] = directories.at(currentChannel).make<TH1D> ("muonD0",channelLabel+" channel: Muon d_{0}; d_{0} [cm] ", 1000, -1, 1);
55    oneDHists_.at(currentChannel)["muonAbsD0"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0",channelLabel+" channel: Muon d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
56    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);
57    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);
58    oneDHists_.at(currentChannel)["muonIso"] = directories.at(currentChannel).make<TH1D> ("muonIso",channelLabel+" channel: Muon Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
59    oneDHists_.at(currentChannel)["numMuons"] = directories.at(currentChannel).make<TH1D> ("numMuons",channelLabel+" channel: Number of Selected Muons; # muons", 10, 0, 10);
60
61
62    //electron histograms
63    oneDHists_.at(currentChannel)["electronPt"]  = directories.at(currentChannel).make<TH1D> ("electronPt",channelLabel+" channel: Electron Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
64    oneDHists_.at(currentChannel)["electronEta"] = directories.at(currentChannel).make<TH1D> ("electronEta",channelLabel+" channel: Electron Eta; #eta", 50, -5, 5);
65    oneDHists_.at(currentChannel)["electronD0"] = directories.at(currentChannel).make<TH1D> ("electronD0",channelLabel+" channel: Electron d_{0}; d_{0} [cm] ", 1000, -1, 1);
66    oneDHists_.at(currentChannel)["electronAbsD0"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0",channelLabel+" channel: Electron d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
67    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);
68    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);
69    oneDHists_.at(currentChannel)["electronIso"] = directories.at(currentChannel).make<TH1D> ("electronIso",channelLabel+" channel: Electron Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
70    oneDHists_.at(currentChannel)["electronFbrem"] = directories.at(currentChannel).make<TH1D> ("electronFbrem",channelLabel+" channel: Electron Brem. Energy Fraction; fbrem", 1000, 0, 2);
71    oneDHists_.at(currentChannel)["numElectrons"] = directories.at(currentChannel).make<TH1D> ("numElectrons",channelLabel+" channel: Number of Selected Electrons; # electrons", 10, 0, 10);
179  
180 +    //book all histograms included in the configuration
181 +    for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
182 +      histogram currentHistogram = histograms.at(currentHistogramIndex);
183 +      int numBinsElements = currentHistogram.bins.size();
184 +      int numInputVariables = currentHistogram.inputVariables.size();
185 +
186 +      if(numBinsElements != 3 && numBinsElements !=6) {
187 +        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
188 +        exit(0);
189 +      }
190 +      else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
191 +        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
192 +        exit(0);
193 +      }
194 +      else if(numBinsElements == 3){
195 +        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));
196 +      }
197 +      else if(numBinsElements == 6){
198 +        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));
199 +
200 +      }
201 +
202 +
203 +      if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
204 +
205 + // bin      particle type
206 + // ---      -------------
207 + //  0        unmatched
208 + //  1        u
209 + //  2        d
210 + //  3        s
211 + //  4        c
212 + //  5        b
213 + //  6        t
214 + //  7        e
215 + //  8        mu
216 + //  9        tau
217 + // 10        nu
218 + // 11        g
219 + // 12        gamma
220 + // 13        Z
221 + // 14        W
222 + // 15        light meson
223 + // 16        K meson
224 + // 17        D meson
225 + // 18        B meson
226 + // 19        light baryon
227 + // 20        strange baryon
228 + // 21        charm baryon
229 + // 22        bottom baryon
230 + // 23        other
231 +
232 +      vector<TString> labelArray;
233 +      labelArray.push_back("unmatched");
234 +      labelArray.push_back("u");
235 +      labelArray.push_back("d");
236 +      labelArray.push_back("s");
237 +      labelArray.push_back("c");
238 +      labelArray.push_back("b");
239 +      labelArray.push_back("t");
240 +      labelArray.push_back("e");
241 +      labelArray.push_back("#mu");
242 +      labelArray.push_back("#tau");
243 +      labelArray.push_back("#nu");
244 +      labelArray.push_back("g");
245 +      labelArray.push_back("#gamma");
246 +      labelArray.push_back("Z");
247 +      labelArray.push_back("W");
248 +      labelArray.push_back("light meson");
249 +      labelArray.push_back("K meson");
250 +      labelArray.push_back("D meson");
251 +      labelArray.push_back("B meson");
252 +      labelArray.push_back("light baryon");
253 +      labelArray.push_back("strange baryon");
254 +      labelArray.push_back("charm baryon");
255 +      labelArray.push_back("bottom baryon");
256 +      labelArray.push_back("other");
257 +
258 +      for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
259 +        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
260 +      }
261 +    }
262 +
263 +    //book a histogram for the number of each object type to be plotted
264 +
265 +    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
266 +      string currentObject = objectsToPlot.at(currentObjectIndex);
267 +      int maxNum = 10;
268 +      if(currentObject == "mcparticles") maxNum = 50;
269 +      else if(currentObject == "primaryvertexs") maxNum = 50;
270 +      else if(currentObject == "muon-muon pairs")
271 +        {
272 +          currentObject = "dimuonPairs";
273 +        }
274 +      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
275 +      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
276 +
277 +      currentObject.at(0) = toupper(currentObject.at(0));
278 +      string histoName = "num" + currentObject;
279 +
280 +      if(histoName == "numPrimaryvertexs"){
281 +        string newHistoName = histoName + "BeforePileupCorrection";
282 +        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);    
283 +        newHistoName = histoName + "AfterPileupCorrection";
284 +        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
285 +      }
286 +      else
287 +        oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
288 +    }
289 +
290 +
291 +    
292  
293  
294      //get list of cuts for this channel
295      vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
296  
78    vector<string> tempInputCollections;
79
80
297      //loop over and parse all cuts
298      for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
83
299        cut tempCut;
300 <     //store input collection for cut
301 <      string inputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
302 <      tempCut.inputCollection = inputCollection;
303 <      
304 <      tempInputCollections.push_back(inputCollection);
305 <      allNecessaryObjects.push_back(inputCollection);
300 >      //store input collection for cut
301 >      string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
302 >      tempCut.inputCollection = tempInputCollection;
303 >      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
304 >        objectsToGet.push_back(tempInputCollection);
305 >        objectsToCut.push_back(tempInputCollection);
306 >      }
307 >      else{//pair of objects, need to add them both to the things to objectsToGet
308 >        int dashIndex = tempInputCollection.find("-");
309 >        int spaceIndex = tempInputCollection.find(" ");
310 >        int secondWordLength = spaceIndex - dashIndex;
311 >        objectsToGet.push_back(tempInputCollection);
312 >        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
313 >        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
314 >        objectsToCut.push_back(tempInputCollection);
315 >        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
316 >        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
317 >
318 >      }
319 >
320 >
321  
322        //split cut string into parts and store them
323        string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
324        std::vector<string> cutStringVector = splitString(cutString);
325 <      tempCut.variable = cutStringVector.at(0);// variable to cut on
326 <      tempCut.comparativeOperator = cutStringVector.at(1);// comparison to make
327 <      tempCut.cutValue = atof(cutStringVector.at(2).c_str());// threshold value to pass cut
325 >      if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
326 >        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
327 >        exit(0);
328 >      }
329 >      tempCut.numSubcuts = (cutStringVector.size()+1)/4;
330 >      for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
331 >        int indexOffset = 4 * subcutIndex;
332 >        string currentVariableString = cutStringVector.at(indexOffset);
333 >        if(currentVariableString.find("(")==std::string::npos){
334 >          tempCut.functions.push_back("");//no function was specified
335 >          tempCut.variables.push_back(currentVariableString);// variable to cut on
336 >        }
337 >        else{
338 >          tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
339 >          string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
340 >          tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
341 >        }
342 >        tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
343 >        tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
344 >        if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
345 >      }
346  
347        //get number of objects required to pass cut for event to pass
348        string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
349        std::vector<string> numberRequiredVector = splitString(numberRequiredString);
350 +      if(numberRequiredVector.size()!=2){
351 +        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
352 +        exit(0);
353 +      }
354  
355 <      // determine number required if comparison contains "="
104 <      int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
105 <      if(numberRequiredVector.at(0) == ">") numberRequiredInt++;
106 <      else if(numberRequiredVector.at(0) == "<") numberRequiredInt--;
107 <
355 >      int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
356        tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
357        tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
358 <      
358 >
359  
360        string tempCutName;
361        if(cuts_.at(currentCut).exists("alias")){
362 <        tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
362 >        tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
363        }
364        else{
365 <        //construct string for cutflow table
366 <        bool plural = numberRequiredInt != 1;
367 <        string collectionString = plural ? inputCollection : inputCollection.substr(0, inputCollection.size()-1);
368 <        string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
369 <        tempCutName = cutName;
365 >        //construct string for cutflow table
366 >        bool plural = numberRequiredInt != 1;
367 >        string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
368 >        string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
369 >        tempCutName = cutName;
370        }
371        tempCut.name = tempCutName;
372  
# Line 126 | Line 374 | OSUAnalysis::OSUAnalysis (const edm::Par
374        tempChannel.cuts.push_back(tempCut);
375  
376  
129    }//end loop over cuts
377  
378 <    //make unique vector of all objects that this channel cuts on (so we know to set flags for those)
132 <    sort( tempInputCollections.begin(), tempInputCollections.end() );
133 <    tempInputCollections.erase( unique( tempInputCollections.begin(), tempInputCollections.end() ), tempInputCollections.end() );
134 <    tempChannel.inputCollections = tempInputCollections;    
378 >    }//end loop over cuts
379  
380      channels.push_back(tempChannel);
381      tempChannel.cuts.clear();
# Line 139 | Line 383 | OSUAnalysis::OSUAnalysis (const edm::Par
383    }//end loop over channels
384  
385  
386 <  //make unique vector of objects we need to cut on (so we make sure to get them from the event)
387 <  sort( allNecessaryObjects.begin(), allNecessaryObjects.end() );
388 <  allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() );
386 >  //make unique vector of objects we need to get from the event
387 >  sort( objectsToGet.begin(), objectsToGet.end() );
388 >  objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
389 >  //make unique vector of objects we need to cut on
390 >  sort( objectsToCut.begin(), objectsToCut.end() );
391 >  objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
392  
393  
394   }
# Line 158 | Line 405 | OSUAnalysis::~OSUAnalysis ()
405   void
406   OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
407   {
408 +
409    // Retrieve necessary collections from the event.
410 <  edm::Handle<BNjetCollection> jets;
411 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end())  
410 >
411 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
412 >    event.getByLabel (triggers_, triggers);
413 >
414 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
415      event.getByLabel (jets_, jets);
416  
417 <  edm::Handle<BNmuonCollection> muons;
167 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end())  
417 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
418      event.getByLabel (muons_, muons);
419  
420 <  edm::Handle<BNelectronCollection> electrons;
171 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end())  
420 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
421      event.getByLabel (electrons_, electrons);
422  
423 <  edm::Handle<BNeventCollection> events;
175 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end())  
423 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
424      event.getByLabel (events_, events);
425  
426 <  edm::Handle<BNtauCollection> taus;
179 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end())  
426 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
427      event.getByLabel (taus_, taus);
428  
429 <  edm::Handle<BNmetCollection> mets;
183 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end())  
429 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
430      event.getByLabel (mets_, mets);
431  
432 <  edm::Handle<BNtrackCollection> tracks;
187 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end())  
432 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
433      event.getByLabel (tracks_, tracks);
434  
435 <  edm::Handle<BNgenjetCollection> genjets;
191 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end())  
435 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
436      event.getByLabel (genjets_, genjets);
437  
438 <  edm::Handle<BNmcparticleCollection> mcparticles;
439 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end())  
440 <    event.getByLabel (mcparticles_, mcparticles);
438 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
439 >        event.getByLabel (mcparticles_, mcparticles);
440 >    
441  
442 <  edm::Handle<BNprimaryvertexCollection> primaryvertexs;
199 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end())  
442 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
443      event.getByLabel (primaryvertexs_, primaryvertexs);
444  
445 <  edm::Handle<BNbxlumiCollection> bxlumis;
203 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end())  
445 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
446      event.getByLabel (bxlumis_, bxlumis);
447  
448 <  edm::Handle<BNphotonCollection> photons;
207 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end())  
448 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
449      event.getByLabel (photons_, photons);
450  
451 <  edm::Handle<BNsuperclusterCollection> superclusters;
211 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end())  
451 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
452      event.getByLabel (superclusters_, superclusters);
453  
454  
455 +  //get pile-up event weight
456 +  double puScaleFactor = 1.00;
457 +  //  if(datasetType_ != "data"){
458 +  // puScaleFactor = puWeight_->at (events->at (0).numTruePV);
459 +    //    cout << puScaleFactor << endl;
460 +  //}
461 +
462 +
463 +
464    //loop over all channels
465  
466    for(uint currentChannelIndex = 0; currentChannelIndex != channels.size(); currentChannelIndex++){
467      channel currentChannel = channels.at(currentChannelIndex);
468  
469 <    flagMap individualFlags;  
470 <    flagMap cumulativeFlags;  
471 <    counterMap passingCounter;
469 >    flagMap individualFlags;
470 >    flagMap cumulativeFlags;
471 >    counterMap passingCounter;
472 >
473 >    bool triggerDecision = true;
474 >    if(currentChannel.triggers.size() != 0){  //triggers specified
475 >      triggerDecision = evaluateTriggers(currentChannel.triggers,triggers.product());
476 >      cutFlows_.at(currentChannelIndex)->at ("trigger") = triggerDecision;
477 >    }
478  
479 +    //loop over all cuts
480  
481  
482  
227    //loop over all cuts
483      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
484        cut currentCut = currentChannel.cuts.at(currentCutIndex);
485  
486 <      for(uint currentObjectIndex = 0; currentObjectIndex != allNecessaryObjects.size(); currentObjectIndex++){
486 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
487 >        string currentObject = objectsToCut.at(currentObjectIndex);
488 >
489 >        //initialize maps to get ready to set cuts
490 >        individualFlags[currentObject].push_back (vector<bool> ());
491 >        cumulativeFlags[currentObject].push_back (vector<bool> ());
492 >
493 >      }
494 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
495 >        string currentObject = objectsToCut.at(currentObjectIndex);
496 >
497 >        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
498  
499 <        string currentObject = allNecessaryObjects.at(currentObjectIndex);
500 <        individualFlags[currentObject].push_back (vector<bool> ());
501 <        cumulativeFlags[currentObject].push_back (vector<bool> ());
502 <
503 <        if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
504 <        else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
505 <        else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
506 <        else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
507 <        else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
508 <        else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
509 <        else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
510 <        else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
511 <        else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
512 <        else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
513 <        else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
514 <        else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
515 <        else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
499 >        if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
500 >        else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
501 >        else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
502 >        else if(currentObject == "events") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,events.product(),"events");
503 >        else if(currentObject == "taus") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,taus.product(),"taus");
504 >        else if(currentObject == "mets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mets.product(),"mets");
505 >        else if(currentObject == "tracks") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,tracks.product(),"tracks");
506 >        else if(currentObject == "genjets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,genjets.product(),"genjets");
507 >        else if(currentObject == "mcparticles") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,mcparticles.product(),"mcparticles");
508 >        else if(currentObject == "primaryvertexs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,primaryvertexs.product(),"primaryvertexs");
509 >        else if(currentObject == "bxlumis") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,bxlumis.product(),"bxlumis");
510 >        else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
511 >        else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
512 >
513 >        
514 >        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
515 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
516 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
517 >                                                                   "muon-muon pairs");
518 >        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
519 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
520 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
521 >                                                                           "electron-electron pairs");
522 >        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
523 >                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
524 >                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
525 >                                                                       "electron-muon pairs");
526  
527  
528        }
# Line 254 | Line 530 | OSUAnalysis::analyze (const edm::Event &
530  
531  
532      }//end loop over all cuts
257    
533  
534  
535      //use cumulative flags to apply cuts at event level
536  
537      bool eventPassedAllCuts = true;
538  
539 +    //apply trigger (true if none were specified)
540 +    eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
541 +
542 +
543 +
544      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
545  
546        //loop over all objects and count how many passed the cumulative selection up to this point
547        cut currentCut = currentChannel.cuts.at(currentCutIndex);
548        int numberPassing = 0;
269      
270      for (uint object = 0; object != cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).size() ; object++)
271          if(cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).at(object)) numberPassing++;
549  
550 +      for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
551 +          if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
552 +      }
553  
554        bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
555        cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
# Line 278 | Line 558 | OSUAnalysis::analyze (const edm::Event &
558  
559      }
560  
561 <    cutFlows_.at(currentChannelIndex)->fillCutFlow();  
282 <  
283 <    string lastCutName = currentChannel.cuts.back().name; //get name of final cut
561 >    cutFlows_.at(currentChannelIndex)->fillCutFlow(puScaleFactor);
562  
563  
564  
565 +    if(!eventPassedAllCuts)continue;
566  
567  
568 <    if(!eventPassedAllCuts)continue;
568 >    
569  
570 +    //set position of primary vertex in event, in order to calculate quantities relative to it
571 +    if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
572 +      vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
573 +      for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
574 +        if(!vertexFlags.at(vertexIndex)) continue;
575 +        chosenPrimaryVertex = & primaryvertexs->at(vertexIndex);
576 +        break;
577 +      }
578 +    }
579 +    else {
580 +      chosenPrimaryVertex = & primaryvertexs->at(0);
581 +    }
582  
583  
293    vector<bool> electronFlags = cumulativeFlags["electrons"].back();
294    int electronCounter = 0;
584  
585 <    for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
586 <      if(!electronFlags.at(electronIndex)) continue;
585 >    //filling histograms
586 >    for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
587 >      histogram currentHistogram = histograms.at(histogramIndex);
588 >
589 >      if(currentHistogram.inputVariables.size() == 1){
590 >        TH1D* histo;
591 >        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
592  
593  
300      electronCounter++;
301      oneDHists_.at(currentChannelIndex)["electronPt"]->Fill (electrons->at(electronIndex).pt);
302      oneDHists_.at(currentChannelIndex)["electronEta"]->Fill (electrons->at(electronIndex).eta);
303      oneDHists_.at(currentChannelIndex)["electronD0"]->Fill (electrons->at(electronIndex).correctedD0Vertex);
304      oneDHists_.at(currentChannelIndex)["electronAbsD0"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex));
305      oneDHists_.at(currentChannelIndex)["electronD0Sig"]->Fill (electrons->at(electronIndex).correctedD0Vertex / electrons->at(electronIndex).tkD0err);
306      oneDHists_.at(currentChannelIndex)["electronAbsD0Sig"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex) / electrons->at(electronIndex).tkD0err);
307      oneDHists_.at(currentChannelIndex)["electronIso"]->Fill (((electrons->at(electronIndex).trackIso + electrons->at(electronIndex).caloIso) / electrons->at(electronIndex).pt));
308      oneDHists_.at(currentChannelIndex)["electronFbrem"]->Fill (electrons->at(electronIndex).fbrem);
309    }
310    oneDHists_.at(currentChannelIndex)["numElectrons"]->Fill (electronCounter);
594  
595 <    
596 <    vector<bool> muonFlags = cumulativeFlags["muons"].back();
597 <    int muonCounter = 0;
595 >        if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
596 >        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
597 >        else if(currentHistogram.inputCollection == "muon-muon pairs")
598 >          {
599 >            fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
600 >                            cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
601 >                            
602 >                            cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
603 >
604 >          }
605 >        else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),puScaleFactor);
606 >        else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
607 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
608 >                                                                                               cumulativeFlags.at("electron-electron pairs").back(),puScaleFactor);
609 >        else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
610 >                                                                                              cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
611 >                                                                                              cumulativeFlags.at("electron-muon pairs").back(),puScaleFactor);
612 >        else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),puScaleFactor);
613 >        else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),puScaleFactor);
614 >        else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),puScaleFactor);
615 >        else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),puScaleFactor);
616 >        else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),puScaleFactor);
617 >        else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),puScaleFactor);
618 >        else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),puScaleFactor);
619 >        else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),puScaleFactor);
620 >        else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),puScaleFactor);
621 >        else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),puScaleFactor);
622 >      }
623 >      else if(currentHistogram.inputVariables.size() == 2){
624 >        TH2D* histo;
625 >        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
626 >
627  
628  
629 <    for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
630 <      if(!muonFlags.at(muonIndex)) continue;
631 <      muonCounter++;
632 <
633 <      oneDHists_.at(currentChannelIndex)["muonPt"]->Fill (muons->at(muonIndex).pt);
634 <      oneDHists_.at(currentChannelIndex)["muonEta"]->Fill (muons->at(muonIndex).eta);
635 <      oneDHists_.at(currentChannelIndex)["muonD0"]->Fill (muons->at(muonIndex).correctedD0Vertex);
636 <      oneDHists_.at(currentChannelIndex)["muonAbsD0"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex));
637 <      oneDHists_.at(currentChannelIndex)["muonD0Sig"]->Fill (muons->at(muonIndex).correctedD0Vertex / muons->at(muonIndex).tkD0err);
638 <      oneDHists_.at(currentChannelIndex)["muonAbsD0Sig"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex) / muons->at(muonIndex).tkD0err);
639 <      oneDHists_.at(currentChannelIndex)["muonIso"]->Fill (((muons->at(muonIndex).trackIso + muons->at(muonIndex).caloIso) / muons->at(muonIndex).pt));
629 >        if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
630 >        else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
631 >        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
632 >                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
633 >                                                                                       cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
634 >        else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),puScaleFactor);
635 >        else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
636 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
637 >                                                                                               cumulativeFlags.at("electron-electron pairs").back(),puScaleFactor);
638 >        else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
639 >                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
640 >                                                                                               cumulativeFlags.at("electron-muon pairs").back(),puScaleFactor);
641 >        else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),puScaleFactor);
642 >        else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),puScaleFactor);
643 >        else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),puScaleFactor);
644 >        else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),puScaleFactor);
645 >        else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),puScaleFactor);
646 >        else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),puScaleFactor);
647 >        else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),puScaleFactor);
648 >        else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),puScaleFactor);
649 >        else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),puScaleFactor);
650 >        else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),puScaleFactor);
651 >      }
652      }
329    oneDHists_.at(currentChannelIndex)["numMuons"]->Fill (muonCounter);
653  
654 +    //fills histograms with the sizes of collections
655 +    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
656 +      string currentObject = objectsToPlot.at(currentObjectIndex);
657 +
658 +      string objectToPlot = "";
659 +
660 +      if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
661 +      else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
662 +      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
663 +      else objectToPlot = currentObject;
664 +      string tempCurrentObject = objectToPlot;
665 +      tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
666 +      string histoName = "num" + tempCurrentObject;
667 +
668 +
669 +
670 +
671 +      //set position of primary vertex in event, in order to calculate quantities relative to it
672 +      if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
673 +        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
674 +        int numToPlot = 0;
675 +        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
676 +          if(lastCutFlags.at(currentFlag)) numToPlot++;
677 +        }
678 +        if(objectToPlot == "primaryvertexs"){
679 +          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
680 +          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
681 +        }
682 +        else
683 +          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,puScaleFactor);        
684 +      }
685 +      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
686 +      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
687 +      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
688 +      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
689 +      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
690 +      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
691 +      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
692 +      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
693 +      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
694 +      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
695 +      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
696 +      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
697 +      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor);
698 +      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor);
699 +      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor);
700 +      else if(objectToPlot == "primaryvertexs"){
701 +        oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
702 +        oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
703 +      }
704 +
705 +    }
706  
707  
708  
709  
710    } //end loop over channel
711  
712 <  masterCutFlow_->fillCutFlow();    
712 >  masterCutFlow_->fillCutFlow(puScaleFactor);
713 >
714  
715  
716   }
717  
718  
719 < bool
719 > bool
720   OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
721  
722 +
723    if(comparison == ">")       return testValue >  cutValue;
724    else if(comparison == ">=") return testValue >= cutValue;
725    else if(comparison == "<")  return testValue <  cutValue;
726    else if(comparison == "<=") return testValue <= cutValue;
727 +  else if(comparison == "==") return testValue == cutValue;
728 +  else if(comparison == "=") return testValue == cutValue;
729 +  else if(comparison == "!=") return testValue != cutValue;
730    else {std::cout << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
731  
732   }
733  
734 + bool
735 + OSUAnalysis::evaluateTriggers (vector<string> triggersToTest, const BNtriggerCollection* triggerCollection){
736 +
737 +  bool triggerDecision = false;
738 +
739 +  for (BNtriggerCollection::const_iterator trigger = triggerCollection->begin (); trigger != triggerCollection->end (); trigger++)
740 +    {
741 +      if(trigger->pass != 1) continue;
742 +      for(uint triggerName = 0; triggerName != triggersToTest.size(); triggerName++)
743 +        {
744 +          if(trigger->name.find(triggersToTest.at(triggerName))!=std::string::npos){
745 +            triggerDecision = true;
746 +          }
747 +        }
748 +    }
749 +  return triggerDecision;
750 +
751 + }
752 +
753   std::vector<std::string>
754   OSUAnalysis::splitString (string inputString){
755  
# Line 363 | Line 762 | OSUAnalysis::splitString (string inputSt
762   }
763  
764   double
765 < OSUAnalysis::valueLookup (const BNjet* object, string variable){
765 > OSUAnalysis::valueLookup (const BNjet* object, string variable, string function){
766  
767    double value = 0.0;
768    if(variable == "energy") value = object->energy;
# Line 481 | Line 880 | OSUAnalysis::valueLookup (const BNjet* o
880    else if(variable == "puJetId_loose_full") value = object->puJetId_loose_full;
881    else if(variable == "puJetId_loose_simple") value = object->puJetId_loose_simple;
882    else if(variable == "puJetId_loose_cutbased") value = object->puJetId_loose_cutbased;
883 <  
884 <  
883 >
884 >
885    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
886 <  
886 >
887 >  value = applyFunction(function, value);
888  
889    return value;
890   }
# Line 492 | Line 892 | OSUAnalysis::valueLookup (const BNjet* o
892  
893  
894   double
895 < OSUAnalysis::valueLookup (const BNmuon* object, string variable){
895 > OSUAnalysis::valueLookup (const BNmuon* object, string variable, string function){
896  
897    double value = 0.0;
898    if(variable == "energy") value = object->energy;
# Line 508 | Line 908 | OSUAnalysis::valueLookup (const BNmuon*
908    else if(variable == "ecalIso") value = object->ecalIso;
909    else if(variable == "hcalIso") value = object->hcalIso;
910    else if(variable == "caloIso") value = object->caloIso;
911 <  else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
911 >  else if(variable == "trackIsDR03") value = object->trackIsoDR03;
912    else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
913    else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
914    else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
# Line 646 | Line 1046 | OSUAnalysis::valueLookup (const BNmuon*
1046    else if(variable == "numberOfMatchedStations") value = object->numberOfMatchedStations;
1047    else if(variable == "time_ndof") value = object->time_ndof;
1048  
1049 +  //user-defined variables
1050 +  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1051 +  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1052 +  else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1053 +  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1054 +  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1055 +
1056 +
1057 +
1058 +  else if(variable == "correctedD0VertexInEBPlus"){
1059 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1060 +    else value = -999;
1061 +  }
1062 +  else if(variable == "correctedD0VertexOutEBPlus"){
1063 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1064 +    else value = -999;
1065 +  }
1066 +  else if(variable == "correctedD0VertexEEPlus"){
1067 +    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex;
1068 +    else value = -999;
1069 +  }
1070 +
1071 +  else if(variable == "correctedD0BeamspotInEBPlus"){
1072 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1073 +    else value = -999;
1074 +  }
1075 +  else if(variable == "correctedD0BeamspotOutEBPlus"){
1076 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1077 +    else value = -999;
1078 +  }
1079 +  else if(variable == "correctedD0BeamspotEEPlus"){
1080 +    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0;
1081 +    else value = -999;
1082 +  }
1083 +
1084 +  else if(variable == "correctedD0VertexInEBMinus"){
1085 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1086 +    else value = -999;
1087 +  }
1088 +  else if(variable == "correctedD0VertexOutEBMinus"){
1089 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1090 +    else value = -999;
1091 +  }
1092 +  else if(variable == "correctedD0VertexEEMinus"){
1093 +    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex;
1094 +    else value = -999;
1095 +  }
1096 +
1097 +  else if(variable == "correctedD0BeamspotInEBMinus"){
1098 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1099 +    else value = -999;
1100 +  }
1101 +  else if(variable == "correctedD0BeamspotOutEBMinus"){
1102 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1103 +    else value = -999;
1104 +  }
1105 +  else if(variable == "correctedD0BeamspotEEMinus"){
1106 +    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0;
1107 +    else value = -999;
1108 +  }
1109 +
1110 +
1111 +  else if(variable == "correctedD0VertexInEBPositiveCharge"){
1112 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1113 +    else value = -999;
1114 +  }
1115 +  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1116 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1117 +    else value = -999;
1118 +  }
1119 +  else if(variable == "correctedD0VertexEEPositiveCharge"){
1120 +    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex;
1121 +    else value = -999;
1122 +  }
1123 +
1124 +  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1125 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1126 +    else value = -999;
1127 +  }
1128 +  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1129 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1130 +    else value = -999;
1131 +  }
1132 +  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1133 +    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0;
1134 +    else value = -999;
1135 +  }
1136 +
1137 +  else if(variable == "correctedD0VertexInEBNegativeCharge"){
1138 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1139 +    else value = -999;
1140 +  }
1141 +  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1142 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1143 +    else value = -999;
1144 +  }
1145 +  else if(variable == "correctedD0VertexEENegativeCharge"){
1146 +    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex;
1147 +    else value = -999;
1148 +  }
1149 +
1150 +  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1151 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1152 +    else value = -999;
1153 +  }
1154 +  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1155 +    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1156 +    else value = -999;
1157 +  }
1158 +  else if(variable == "correctedD0BeamspotEENegativeCharge"){
1159 +    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0;
1160 +    else value = -999;
1161 +  }
1162 +
1163 +
1164 +  else if(variable == "tightID") {
1165 +    value = object->isGlobalMuon > 0            \
1166 +      && object->isPFMuon > 0                   \
1167 +      && object->normalizedChi2 < 10            \
1168 +      && object->numberOfValidMuonHits > 0      \
1169 +      && object->numberOfMatchedStations > 1    \
1170 +      && fabs(object->correctedD0Vertex) < 0.2  \
1171 +      && fabs(object->correctedDZ) < 0.5        \
1172 +      && object->numberOfValidPixelHits > 0             \
1173 +      && object->numberOfLayersWithMeasurement > 5;
1174 +  }
1175 +  else if(variable == "tightIDdisplaced"){
1176 +    value = object->isGlobalMuon > 0            \
1177 +      && object->normalizedChi2 < 10            \
1178 +      && object->numberOfValidMuonHits > 0      \
1179 +      && object->numberOfMatchedStations > 1    \
1180 +      && object->numberOfValidPixelHits > 0     \
1181 +      && object->numberOfLayersWithMeasurement > 5;
1182 +  }
1183 +
1184 +
1185 +
1186 +  else if(variable == "genMatchedId"){
1187 +    int index = getGenMatchedParticleIndex(object);
1188 +    if(index == -1) value = 0;
1189 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1190 +  }
1191 +  else if(variable == "genMatchedMotherId"){
1192 +    int index = getGenMatchedParticleIndex(object);
1193 +    if(index == -1) value = 0;
1194 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1195 +  }
1196 +  else if(variable == "genMatchedGrandmotherId"){
1197 +    int index = getGenMatchedParticleIndex(object);
1198 +    if(index == -1) value = 0;
1199 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1200 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1201 +      if(motherIndex == -1) value = 0;
1202 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1203 +    }
1204 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1205 +  }
1206 +
1207 +
1208 +
1209    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1210 <  
1210 >
1211 >  value = applyFunction(function, value);
1212  
1213    return value;
1214   }
1215  
1216  
1217   double
1218 < OSUAnalysis::valueLookup (const BNelectron* object, string variable){
1218 > OSUAnalysis::valueLookup (const BNelectron* object, string variable, string function){
1219  
1220    double value = 0.0;
1221    if(variable == "energy") value = object->energy;
# Line 814 | Line 1375 | OSUAnalysis::valueLookup (const BNelectr
1375    else if(variable == "eidHyperTight4MC") value = object->eidHyperTight4MC;
1376    else if(variable == "passConvVeto") value = object->passConvVeto;
1377  
1378 +  //user-defined variables
1379 +  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1380 +  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1381 +  else if(variable == "detIso") value = (object->trackIso) / object->pt;
1382 +  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1383 +
1384 +
1385 +  else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
1386 +    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
1387 +    else value = -999;
1388 +  }
1389 +  else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){
1390 +    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex;
1391 +    else value = -999;
1392 +  }
1393 +
1394 +
1395 +  else if(variable == "correctedD0VertexInEBPlus"){
1396 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1397 +    else value = -999;
1398 +  }
1399 +  else if(variable == "correctedD0VertexOutEBPlus"){
1400 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
1401 +    else value = -999;
1402 +  }
1403 +  else if(variable == "correctedD0VertexEEPlus"){
1404 +    if(object->isEE && object->eta > 0) value = object->correctedD0Vertex;
1405 +    else value = -999;
1406 +  }
1407 +
1408 +  else if(variable == "correctedD0BeamspotInEBPlus"){
1409 +    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
1410 +    else value = -999;
1411 +  }
1412 +  else if(variable == "correctedD0BeamspotOutEBPlus"){
1413 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
1414 +    else value = -999;
1415 +  }
1416 +  else if(variable == "correctedD0BeamspotEEPlus"){
1417 +    if(object->isEE && object->eta > 0) value = object->correctedD0;
1418 +    else value = -999;
1419 +  }
1420 +
1421 +  else if(variable == "correctedD0VertexInEBMinus"){
1422 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1423 +    else value = -999;
1424 +  }
1425 +  else if(variable == "correctedD0VertexOutEBMinus"){
1426 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
1427 +    else value = -999;
1428 +  }
1429 +  else if(variable == "correctedD0VertexEEMinus"){
1430 +    if(object->isEE && object->eta < 0) value = object->correctedD0Vertex;
1431 +    else value = -999;
1432 +  }
1433 +
1434 +  else if(variable == "correctedD0BeamspotInEBMinus"){
1435 +    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
1436 +    else value = -999;
1437 +  }
1438 +  else if(variable == "correctedD0BeamspotOutEBMinus"){
1439 +    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
1440 +    else value = -999;
1441 +  }
1442 +  else if(variable == "correctedD0BeamspotEEMinus"){
1443 +    if(object->isEE && object->eta < 0) value = object->correctedD0;
1444 +    else value = -999;
1445 +  }
1446 +
1447 +
1448 +  else if(variable == "correctedD0VertexInEBPositiveCharge"){
1449 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1450 +    else value = -999;
1451 +  }
1452 +  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
1453 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
1454 +    else value = -999;
1455 +  }
1456 +  else if(variable == "correctedD0VertexEEPositiveCharge"){
1457 +    if(object->isEE && object->charge > 0) value = object->correctedD0Vertex;
1458 +    else value = -999;
1459 +  }
1460 +
1461 +  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
1462 +    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
1463 +    else value = -999;
1464 +  }
1465 +  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
1466 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
1467 +    else value = -999;
1468 +  }
1469 +  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
1470 +    if(object->isEE && object->charge > 0) value = object->correctedD0;
1471 +    else value = -999;
1472 +  }
1473 +
1474 +  else if(variable == "correctedD0VertexInEBNegativeCharge"){
1475 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1476 +    else value = -999;
1477 +  }
1478 +  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
1479 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
1480 +    else value = -999;
1481 +  }
1482 +  else if(variable == "correctedD0VertexEENegativeCharge"){
1483 +    if(object->isEE && object->charge < 0) value = object->correctedD0Vertex;
1484 +    else value = -999;
1485 +  }
1486 +
1487 +  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
1488 +    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
1489 +    else value = -999;
1490 +  }
1491 +  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
1492 +    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
1493 +    else value = -999;
1494 +  }
1495 +  else if(variable == "correctedD0BeamspotEENegativeCharge"){
1496 +    if(object->isEE && object->charge < 0) value = object->correctedD0;
1497 +    else value = -999;
1498 +  }
1499 +
1500 +
1501 +  else if(variable == "tightIDdisplaced"){
1502 +    if (object->isEB)
1503 +      {
1504 +        value = fabs(object->delEtaIn) < 0.004 \
1505 +          && fabs (object->delPhiIn) < 0.03 \
1506 +          && object->scSigmaIEtaIEta < 0.01 \
1507 +          && object->hadOverEm < 0.12 \
1508 +          && object->absInvEMinusInvPin < 0.05;
1509 +      }
1510 +    else
1511 +      {
1512 +        value = fabs (object->delEtaIn) < 0.005 \
1513 +          && fabs (object->delPhiIn) < 0.02 \
1514 +          && object->scSigmaIEtaIEta < 0.03 \
1515 +          && object->hadOverEm < 0.10 \
1516 +          && object->absInvEMinusInvPin < 0.05;
1517 +      }
1518 +  }
1519 +
1520 +  else if(variable == "genMatchedId"){
1521 +    int index = getGenMatchedParticleIndex(object);
1522 +    if(index == -1) value = 0;
1523 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1524 +  }
1525 +  else if(variable == "genMatchedMotherId"){
1526 +    int index = getGenMatchedParticleIndex(object);
1527 +    if(index == -1) value = 0;
1528 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1529 +  }
1530 +  else if(variable == "genMatchedGrandmotherId"){
1531 +    int index = getGenMatchedParticleIndex(object);
1532 +    if(index == -1) value = 0;
1533 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1534 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1535 +      if(motherIndex == -1) value = 0;
1536 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1537 +    }
1538 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1539 +  }
1540 +
1541 +
1542  
1543    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1544 <  
1544 >
1545 >  value = applyFunction(function, value);
1546  
1547    return value;
1548   }
1549  
1550  
1551   double
1552 < OSUAnalysis::valueLookup (const BNevent* object, string variable){
1552 > OSUAnalysis::valueLookup (const BNevent* object, string variable, string function){
1553  
1554    double value = 0.0;
1555  
# Line 893 | Line 1619 | OSUAnalysis::valueLookup (const BNevent*
1619    else if(variable == "evt") value = object->evt;
1620  
1621    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1622 <  
1622 >
1623 >  value = applyFunction(function, value);
1624  
1625    return value;
1626   }
1627  
1628   double
1629 < OSUAnalysis::valueLookup (const BNtau* object, string variable){
1629 > OSUAnalysis::valueLookup (const BNtau* object, string variable, string function){
1630  
1631    double value = 0.0;
1632  
# Line 944 | Line 1671 | OSUAnalysis::valueLookup (const BNtau* o
1671    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1672  
1673  
1674 +  else if(variable == "genMatchedId"){
1675 +    int index = getGenMatchedParticleIndex(object);
1676 +    if(index == -1) value = 0;
1677 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1678 +  }
1679 +  else if(variable == "genMatchedMotherId"){
1680 +    int index = getGenMatchedParticleIndex(object);
1681 +    if(index == -1) value = 0;
1682 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1683 +  }
1684 +  else if(variable == "genMatchedGrandmotherId"){
1685 +    int index = getGenMatchedParticleIndex(object);
1686 +    if(index == -1) value = 0;
1687 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1688 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1689 +      if(motherIndex == -1) value = 0;
1690 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1691 +    }
1692 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1693 +  }
1694 +
1695 +
1696    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1697 <  
1697 >
1698 >  value = applyFunction(function, value);
1699  
1700    return value;
1701   }
1702  
1703   double
1704 < OSUAnalysis::valueLookup (const BNmet* object, string variable){
1704 > OSUAnalysis::valueLookup (const BNmet* object, string variable, string function){
1705  
1706    double value = 0.0;
1707  
# Line 1016 | Line 1766 | OSUAnalysis::valueLookup (const BNmet* o
1766    else if(variable == "pfT1jet10phi") value = object->pfT1jet10phi;
1767  
1768    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1769 <  
1769 >
1770 >  value = applyFunction(function, value);
1771  
1772    return value;
1773   }
1774  
1775   double
1776 < OSUAnalysis::valueLookup (const BNtrack* object, string variable){
1776 > OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){
1777  
1778    double value = 0.0;
1779 <
1779 >  double pMag = sqrt(object->pt * object->pt +
1780 >                         object->pz * object->pz);
1781 >  
1782    if(variable == "pt") value = object->pt;
1783    else if(variable == "px") value = object->px;
1784    else if(variable == "py") value = object->py;
# Line 1045 | Line 1798 | OSUAnalysis::valueLookup (const BNtrack*
1798    else if(variable == "isHighPurity") value = object->isHighPurity;
1799  
1800  
1801 +  //additional BNs info for disappTrks  
1802 +  else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
1803 +  else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
1804 +  else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
1805 +  else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
1806 +  else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
1807 +  else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
1808 +  else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
1809 +  else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
1810 +  else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
1811 +  else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
1812 +  //user defined variables
1813 +  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;
1814 +  else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
1815 +  else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
1816 +  else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
1817 +  else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
1818 +  else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
1819 +  else if(variable == "ptErrorByPt") value = (object->ptError/object->pt);
1820 +  else if(variable == "ptError") value = object->ptError;
1821 +  else if(variable == "ptRes") value = getTrkPtRes(object);  
1822 +
1823 +
1824 +  else if(variable == "genMatchedId"){
1825 +    int index = getGenMatchedParticleIndex(object);
1826 +    if(index == -1) value = 0;
1827 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1828 +  }
1829 +  else if(variable == "genMatchedMotherId"){
1830 +    int index = getGenMatchedParticleIndex(object);
1831 +    if(index == -1) value = 0;
1832 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1833 +  }
1834 +  else if(variable == "genMatchedGrandmotherId"){
1835 +    int index = getGenMatchedParticleIndex(object);
1836 +    if(index == -1) value = 0;
1837 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1838 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1839 +      if(motherIndex == -1) value = 0;
1840 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1841 +    }
1842 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1843 +  }
1844 +
1845 +
1846 +
1847    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1848 <  
1848 >
1849 >  value = applyFunction(function, value);
1850  
1851    return value;
1852   }
1853  
1854   double
1855 < OSUAnalysis::valueLookup (const BNgenjet* object, string variable){
1855 > OSUAnalysis::valueLookup (const BNgenjet* object, string variable, string function){
1856  
1857    double value = 0.0;
1858  
# Line 1073 | Line 1873 | OSUAnalysis::valueLookup (const BNgenjet
1873  
1874  
1875    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1876 <  
1876 >
1877 >  value = applyFunction(function, value);
1878  
1879    return value;
1880   }
1881  
1882   double
1883 < OSUAnalysis::valueLookup (const BNmcparticle* object, string variable){
1883 > OSUAnalysis::valueLookup (const BNmcparticle* object, string variable, string function){
1884  
1885    double value = 0.0;
1886  
# Line 1168 | Line 1969 | OSUAnalysis::valueLookup (const BNmcpart
1969    else if(variable == "grandMother11Status") value = object->grandMother11Status;
1970    else if(variable == "grandMother11Charge") value = object->grandMother11Charge;
1971  
1972 +  //user-defined variables
1973 +  else if (variable == "d0"){
1974 +    double vx = object->vx - chosenPrimaryVertex->x,
1975 +      vy = object->vy - chosenPrimaryVertex->y,
1976 +      px = object->px,
1977 +      py = object->py,
1978 +      pt = object->pt;
1979 +    value = (-vx * py + vy * px) / pt;
1980 +  }
1981 +  else if (variable == "dz"){
1982 +    double vx = object->vx - chosenPrimaryVertex->x,
1983 +      vy = object->vy - chosenPrimaryVertex->y,
1984 +      vz = object->vz - chosenPrimaryVertex->z,
1985 +      px = object->px,
1986 +      py = object->py,
1987 +      pz = object->pz,
1988 +      pt = object->pt;
1989 +    value = vz - (vx * px + vy * py)/pt * (pz/pt);
1990 +  }
1991 +  else if(variable == "v0"){
1992 +    value = sqrt(object->vx*object->vx + object->vy*object->vy);
1993 +  }
1994 +  else if(variable == "deltaV0"){
1995 +    value = sqrt((object->vx-chosenPrimaryVertex->x)*(object->vx-chosenPrimaryVertex->x) + (object->vy-chosenPrimaryVertex->y)*(object->vy-chosenPrimaryVertex->y));
1996 +  }
1997 +  else if (variable == "deltaVx"){
1998 +    value = object->vx - chosenPrimaryVertex->x;
1999 +  }
2000 +  else if (variable == "deltaVy"){
2001 +    value = object->vy - chosenPrimaryVertex->y;
2002 +  }
2003 +  else if (variable == "deltaVz"){
2004 +    value = object->vz - chosenPrimaryVertex->z;
2005 +  }
2006 +
2007 +
2008    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2009 <  
2009 >
2010 >  value = applyFunction(function, value);
2011  
2012    return value;
2013   }
2014  
2015   double
2016 < OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable){
2016 > OSUAnalysis::valueLookup (const BNprimaryvertex* object, string variable, string function){
2017  
2018    double value = 0.0;
2019  
# Line 1195 | Line 2033 | OSUAnalysis::valueLookup (const BNprimar
2033  
2034  
2035    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2036 <  
2036 >
2037 >  value = applyFunction(function, value);
2038  
2039    return value;
2040   }
2041  
2042   double
2043 < OSUAnalysis::valueLookup (const BNbxlumi* object, string variable){
2043 > OSUAnalysis::valueLookup (const BNbxlumi* object, string variable, string function){
2044  
2045    double value = 0.0;
2046  
# Line 1211 | Line 2050 | OSUAnalysis::valueLookup (const BNbxlumi
2050  
2051  
2052    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2053 <  
2053 >
2054 >  value = applyFunction(function, value);
2055  
2056    return value;
2057   }
2058  
2059   double
2060 < OSUAnalysis::valueLookup (const BNphoton* object, string variable){
2060 > OSUAnalysis::valueLookup (const BNphoton* object, string variable, string function){
2061  
2062    double value = 0.0;
2063  
# Line 1292 | Line 2132 | OSUAnalysis::valueLookup (const BNphoton
2132    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
2133  
2134  
2135 +  else if(variable == "genMatchedId"){
2136 +    int index = getGenMatchedParticleIndex(object);
2137 +    if(index == -1) value = 0;
2138 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
2139 +  }
2140 +  else if(variable == "genMatchedMotherId"){
2141 +    int index = getGenMatchedParticleIndex(object);
2142 +    if(index == -1) value = 0;
2143 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
2144 +  }
2145 +  else if(variable == "genMatchedGrandmotherId"){
2146 +    int index = getGenMatchedParticleIndex(object);
2147 +    if(index == -1) value = 0;
2148 +    else if(fabs(mcparticles->at(index).motherId) == 15){
2149 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
2150 +      if(motherIndex == -1) value = 0;
2151 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
2152 +    }
2153 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
2154 +  }
2155 +
2156 +
2157    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2158 <  
2158 >
2159 >  value = applyFunction(function, value);
2160  
2161    return value;
2162   }
2163  
2164   double
2165 < OSUAnalysis::valueLookup (const BNsupercluster* object, string variable){
2165 > OSUAnalysis::valueLookup (const BNsupercluster* object, string variable, string function){
2166  
2167    double value = 0.0;
2168  
# Line 1314 | Line 2177 | OSUAnalysis::valueLookup (const BNsuperc
2177  
2178  
2179    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2180 +
2181 +  value = applyFunction(function, value);
2182 +
2183 +  return value;
2184 + }
2185 +
2186 +
2187 + double
2188 + OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function){
2189 +
2190 +  double value = 0.0;
2191 +
2192 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2193 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2194 +  else if(variable == "invMass"){
2195 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2196 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2197 +
2198 +    value = (fourVector1 + fourVector2).M();}
2199 +
2200 +  else if(variable == "threeDAngle")
2201 +    {
2202 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2203 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2204 +      value = (threeVector1.Angle(threeVector2));
2205 +    }
2206 +  else if(variable == "alpha")
2207 +    {
2208 +      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
2209 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2210 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2211 +      value = (pi-threeVector1.Angle( threeVector2));
2212 +    }
2213 +
2214 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2215 +
2216 +
2217 +  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2218 +  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2219 +  else if(variable == "d0Sign"){
2220 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2221 +    if(d0Sign < 0) value = -0.5;
2222 +    else if (d0Sign > 0) value = 0.5;
2223 +    else value = -999;
2224 +  }
2225 +  else if(variable == "muon1CorrectedD0Vertex"){
2226 +    value = object1->correctedD0Vertex;
2227 +  }
2228 +  else if(variable == "muon2CorrectedD0Vertex"){
2229 +    value = object2->correctedD0Vertex;
2230 +  }
2231 + else if(variable == "muon1timeAtIpInOut"){
2232 +    value = object1->timeAtIpInOut;
2233 +  }
2234 + else if(variable == "muon2timeAtIpInOut"){
2235 +    value = object2->timeAtIpInOut;
2236 +  }
2237 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2238 +
2239 +
2240 +
2241 +
2242 +
2243 +
2244 +
2245 +
2246 +  value = applyFunction(function, value);
2247 +
2248 +  return value;
2249 + }
2250 +
2251 + double
2252 + OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function){
2253 +
2254 +  double value = 0.0;
2255 +
2256 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2257 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2258 +  else if(variable == "invMass"){
2259 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2260 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2261 +
2262 +    value = (fourVector1 + fourVector2).M();}
2263 + else if(variable == "threeDAngle")
2264 +    {
2265 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2266 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2267 +      value = (threeVector1.Angle(threeVector2));
2268 +    }
2269    
2270  
2271 +
2272 +
2273 +  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2274 +  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2275 +
2276 +  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2277 +
2278 +
2279 +
2280 +  else if(variable == "d0Sign"){
2281 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2282 +    if(d0Sign < 0) value = -0.5;
2283 +    else if (d0Sign > 0) value = 0.5;
2284 +    else value = -999;
2285 +  }
2286 +  else if(variable == "electron1CorrectedD0Vertex"){
2287 +    value = object1->correctedD0Vertex;
2288 +  }
2289 +  else if(variable == "electron2CorrectedD0Vertex"){
2290 +    value = object2->correctedD0Vertex;
2291 +  }
2292 +
2293 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2294 +
2295 +  value = applyFunction(function, value);
2296 +
2297 +  return value;
2298 + }
2299 +
2300 + double
2301 + OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function){
2302 +
2303 +  double value = 0.0;
2304 +
2305 +  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2306 +  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2307 +  else if(variable == "invMass"){
2308 +    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2309 +    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2310 +
2311 +    value = (fourVector1 + fourVector2).M();}
2312 + else if(variable == "threeDAngle")
2313 +    {
2314 +      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2315 +      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2316 +      value = (threeVector1.Angle(threeVector2));
2317 +    }
2318 +
2319 +
2320 +  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2321 +  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2322 +  else if(variable == "d0Sign"){
2323 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2324 +    if(d0Sign < 0) value = -0.5;
2325 +    else if (d0Sign > 0) value = 0.5;
2326 +    else value = -999;
2327 +  }
2328 +  else if(variable == "electronCorrectedD0Vertex"){
2329 +    value = object1->correctedD0Vertex;
2330 +  }
2331 +  else if(variable == "muonCorrectedD0Vertex"){
2332 +    value = object2->correctedD0Vertex;
2333 +  }
2334 +  else if(variable == "electronCorrectedD0"){
2335 +    value = object1->correctedD0;
2336 +  }
2337 +  else if(variable == "muonCorrectedD0"){
2338 +    value = object2->correctedD0;
2339 +  }
2340 +
2341 +
2342 +  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2343 +
2344 +  value = applyFunction(function, value);
2345 +
2346 +  return value;
2347 + }
2348 +
2349 +
2350 + // Calculate the number of tracks in cone of DeltaR<0.5 around track1.  
2351 + // Return true iff no other tracks are found in this cone.  
2352 + int
2353 + OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
2354 +  for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
2355 +    if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.  
2356 +    double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
2357 +    if(deltaRtrk < 0.5) return 0;
2358 +  }
2359 +  return 1;
2360 +  
2361 + }
2362 +
2363 +
2364 + double
2365 + OSUAnalysis::getTrkPtRes (const BNtrack* track1){
2366 +
2367 +  double ptTrue = getTrkPtTrue(track1, mcparticles.product());
2368 +  double PtRes = (track1->pt - ptTrue) / ptTrue;
2369 +
2370 +  return PtRes;
2371 +  
2372 + }
2373 +
2374 +
2375 + double
2376 + OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
2377 +  double value = -99;
2378 +  double genDeltaRLowest = 999;
2379 +
2380 +  for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
2381 +    double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
2382 +    if (genDeltaRtemp < genDeltaRLowest) {
2383 +      genDeltaRLowest = genDeltaRtemp;
2384 +      if (genDeltaRLowest < 0.05) {   // Only consider it truth-matched if DeltaR<0.15.  
2385 +        double ptTrue = genPart->pt;
2386 +        value = ptTrue;
2387 +      }
2388 +    }
2389 +  }
2390 +
2391    return value;
2392 +
2393   }
2394  
2395 + //creates a map of the dead Ecal channels in the barrel and endcap
2396 + //to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
2397 + //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
2398 + void
2399 + OSUAnalysis::WriteDeadEcal (){
2400 +  double etaEcal, phiEcal;
2401 +  ifstream DeadEcalFile(deadEcalFile_);
2402 +  if(!DeadEcalFile) {
2403 +    cout << "Error: DeadEcalFile has not been found." << endl;
2404 +    return;
2405 +  }
2406 +  if(DeadEcalVec.size()!= 0){
2407 +    cout << "Error: DeadEcalVec has a nonzero size" << endl;
2408 +    return;
2409 +  }
2410 +  while(!DeadEcalFile.eof())
2411 +    {
2412 +      DeadEcalFile >> etaEcal >> phiEcal;
2413 +      DeadEcal newChan;
2414 +      newChan.etaEcal = etaEcal;
2415 +      newChan.phiEcal = phiEcal;
2416 +      DeadEcalVec.push_back(newChan);
2417 +    }
2418 +  if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl;
2419 + }
2420  
2421 < template <class InputCollection>
2421 > //if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
2422 > int
2423 > OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
2424 >  double deltaRLowest = 999;
2425 >  int value = 0;
2426 >  if (DeadEcalVec.size() == 0) WriteDeadEcal();
2427 >  for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
2428 >    double eta = ecal->etaEcal;
2429 >    double phi = ecal->phiEcal;
2430 >    double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
2431 >    if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
2432 >  }
2433 >  if (deltaRLowest<0.05) {value = 1;}
2434 >  else {value = 0;}
2435 >  return value;
2436 > }
2437 >
2438 >
2439 >
2440 >
2441 > double
2442 > OSUAnalysis::applyFunction(string function, double value){
2443 >
2444 >  if(function == "abs") value = fabs(value);
2445 >  else if(function == "fabs") value = fabs(value);
2446 >  else if(function == "log10") value = log10(value);
2447 >  else if(function == "log") value = log10(value);
2448 >
2449 >  else if(function == "") value = value;
2450 >  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
2451 >
2452 >  return value;
2453 >
2454 > }
2455 >
2456 >
2457 > template <class InputCollection>
2458   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2459  
2460  
2461    for (uint object = 0; object != inputCollection->size(); object++){
2462 <    
2462 >
2463 >
2464      bool decision = true;//object passes if this cut doesn't cut on that type of object
1330    
1331    if(currentCut.inputCollection == inputType){
1332      
1333      double value = valueLookup(&inputCollection->at(object), currentCut.variable);
2465  
2466 <      decision = evaluateComparison(value,currentCut.comparativeOperator,currentCut.cutValue);
2466 >
2467 >    if(currentCut.inputCollection == inputType){
2468 >
2469 >      vector<bool> subcutDecisions;
2470 >      for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2471 >        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2472 >        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2473 >      }
2474 >      if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2475 >      else{
2476 >        bool tempDecision = true;
2477 >        for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2478 >          if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2479 >            tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2480 >          else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2481 >            tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
2482 >        }
2483 >        decision = tempDecision;
2484 >      }
2485      }
2486 <    individualFlags[inputType].at(currentCutIndex).push_back(decision);
2487 <    
2486 >    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2487 >
2488  
2489      //set flags for objects that pass each cut AND all the previous cuts
2490      bool previousCumulativeFlag = true;
2491      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2492 <      if(previousCumulativeFlag && individualFlags[inputType].at(previousCutIndex).at(object)) previousCumulativeFlag = true;
2492 >      if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(object)) previousCumulativeFlag = true;
2493        else{ previousCumulativeFlag = false; break;}
2494      }
2495 <    cumulativeFlags[inputType].at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2495 >    cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2496 >
2497 >
2498 >  }
2499 >
2500 > }
2501 >
2502 >
2503 > template <class InputCollection1, class InputCollection2>
2504 > void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
2505 >                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
2506 >
2507 >
2508 >  bool sameObjects = false;
2509 >  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2510 >
2511 >
2512 >  int counter = 0;  
2513 >  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2514 >    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2515 >      
2516 >      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2517 >
2518 >
2519 >      bool decision = true;//object passes if this cut doesn't cut on that type of object
2520 >
2521 >      if(currentCut.inputCollection == inputType){
2522 >
2523 >        vector<bool> subcutDecisions;
2524 >        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
2525 >          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
2526 >          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
2527 >        }
2528 >
2529 >        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
2530 >        else{
2531 >          bool tempDecision = true;
2532 >          for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
2533 >            if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
2534 >              tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
2535 >            else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
2536 >              tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
2537 >          }
2538 >          decision = tempDecision;
2539 >        }
2540 >      }
2541 >      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2542 >      
2543 >      //set flags for objects that pass each cut AND all the previous cuts
2544 >      bool previousCumulativeFlag = true;
2545 >      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
2546 >        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
2547 >        else{ previousCumulativeFlag = false; break;}
2548 >      }
2549 >      //apply flags for the components of the composite object as well
2550 >      bool currentCumulativeFlag = true;
2551 >      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
2552 >      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
2553 >      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
2554 >      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
2555 >      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
2556 >
2557 >      counter++;      
2558 >    }
2559      
2560    }
2561 <  
2561 >
2562 >
2563 > }
2564 >
2565 >
2566 > template <class InputCollection>
2567 > void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double puScaleFactor){
2568 >
2569 >  for (uint object = 0; object != inputCollection->size(); object++){
2570 >
2571 >    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2572 >
2573 >    string currentString = parameters.inputVariables.at(0);
2574 >    string inputVariable = "";
2575 >    string function = "";
2576 >    if(currentString.find("(")==std::string::npos){
2577 >      inputVariable = currentString;// variable to cut on                                                                                                                                        
2578 >    }
2579 >    else{
2580 >      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2581 >      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2582 >      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2583 >    }
2584 >
2585 >    double value = valueLookup(&inputCollection->at(object), inputVariable, function);
2586 >    histo->Fill(value,puScaleFactor);
2587 >
2588 >  }
2589 >
2590   }
2591  
2592 + template <class InputCollection1, class InputCollection2>
2593 + void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double puScaleFactor){
2594 +
2595 +  bool sameObjects = false;
2596 +  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2597  
2598 +  int pairCounter = 0;  
2599 +  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2600 +    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2601 +      
2602 +      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2603 +
2604 +      //only take objects which have passed all cuts and pairs which have passed all cuts
2605 +      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2606 +      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2607 +      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2608 +
2609 +      string currentString = parameters.inputVariables.at(0);
2610 +      string inputVariable = "";
2611 +      string function = "";
2612 +      if(currentString.find("(")==std::string::npos){
2613 +        inputVariable = currentString;// variable to cut on                                                                                                          
2614 +      }
2615 +      else{
2616 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2617 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2618 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2619 +      }
2620 +      
2621 +      double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2622 +      histo->Fill(value,puScaleFactor);
2623 +
2624 +      pairCounter++;      
2625 +    }
2626 +  }
2627 +
2628 + }
2629 +
2630 +
2631 + template <class InputCollection>
2632 + void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double puScaleFactor){
2633 +
2634 +  for (uint object = 0; object != inputCollection->size(); object++){
2635  
2636 +    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
2637 +
2638 +    string currentString = parameters.inputVariables.at(0);
2639 +    string inputVariable = "";
2640 +    string function = "";
2641 +    if(currentString.find("(")==std::string::npos){
2642 +      inputVariable = currentString;// variable to cut on                                                                                                                                        
2643 +    }
2644 +    else{
2645 +      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2646 +      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2647 +      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2648 +    }
2649 +    double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
2650 +
2651 +    currentString = parameters.inputVariables.at(1);
2652 +    inputVariable = "";
2653 +    function = "";
2654 +    if(currentString.find("(")==std::string::npos){
2655 +      inputVariable = currentString;// variable to cut on                                                                                                                                        
2656 +    }
2657 +    else{
2658 +      function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2659 +      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2660 +      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2661 +    }
2662 +
2663 +    double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
2664 +
2665 +    histo->Fill(valueX,valueY,puScaleFactor);
2666 +
2667 +  }
2668 +
2669 + }
2670 +
2671 + template <class InputCollection1, class InputCollection2>
2672 + void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double puScaleFactor){
2673 +
2674 +  bool sameObjects = false;
2675 +  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2676 +
2677 +  int pairCounter = 0;  
2678 +  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2679 +    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2680 +      
2681 +      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
2682 +
2683 +      //only take objects which have passed all cuts and pairs which have passed all cuts
2684 +      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2685 +      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2686 +      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2687 +
2688 +      string currentString = parameters.inputVariables.at(0);
2689 +      string inputVariable = "";
2690 +      string function = "";
2691 +      if(currentString.find("(")==std::string::npos){
2692 +        inputVariable = currentString;// variable to cut on                                                                                                          
2693 +      }
2694 +      else{
2695 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2696 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2697 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2698 +      }
2699 +      double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2700 +
2701 +      currentString = parameters.inputVariables.at(1);
2702 +      inputVariable = "";
2703 +      function = "";
2704 +      if(currentString.find("(")==std::string::npos){
2705 +        inputVariable = currentString;// variable to cut on                                                                                                          
2706 +      }
2707 +      else{
2708 +        function = currentString.substr(0,currentString.find("("));//function comes before the "("                                                                                          
2709 +        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string                                                                                                    
2710 +        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"                                                                                                    
2711 +      }
2712 +      double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2713 +
2714 +
2715 +      histo->Fill(valueX,valueY,puScaleFactor);
2716 +
2717 +      pairCounter++;      
2718 +
2719 +    }
2720 +  }
2721 +
2722 + }
2723 +
2724 +
2725 + template <class InputObject>
2726 + int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
2727 +
2728 +  int bestMatchIndex = -1;
2729 +  double bestMatchDeltaR = 999;
2730 +
2731 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2732 +
2733 +    double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
2734 +    if(currentDeltaR > 0.05) continue;
2735 +
2736 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2737 +      bestMatchIndex = mcparticle - mcparticles->begin();
2738 +      bestMatchDeltaR = currentDeltaR;
2739 +    }
2740 +
2741 +  }
2742 +
2743 +  return bestMatchIndex;
2744 +
2745 + }
2746 +
2747 +
2748 + int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
2749 +
2750 +  int bestMatchIndex = -1;
2751 +  double bestMatchDeltaR = 999;
2752 +
2753 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2754 +
2755 +    if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
2756 +
2757 +    double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
2758 +    if(currentDeltaR > 0.05) continue;
2759 +
2760 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2761 +      bestMatchIndex = mcparticle - mcparticles->begin();
2762 +      bestMatchDeltaR = currentDeltaR;
2763 +    }
2764 +
2765 +  }
2766 +
2767 +  return bestMatchIndex;
2768 + }
2769 +
2770 +
2771 + // bin      particle type
2772 + // ---      -------------
2773 + //  0        unmatched
2774 + //  1        u
2775 + //  2        d
2776 + //  3        s
2777 + //  4        c
2778 + //  5        b
2779 + //  6        t
2780 + //  7        e
2781 + //  8        mu
2782 + //  9        tau
2783 + // 10        nu
2784 + // 11        g
2785 + // 12        gamma
2786 + // 13        Z
2787 + // 14        W
2788 + // 15        light meson
2789 + // 16        K meson
2790 + // 17        D meson
2791 + // 18        B meson
2792 + // 19        light baryon
2793 + // 20        strange baryon
2794 + // 21        charm baryon
2795 + // 22        bottom baryon
2796 + // 23        other
2797 +
2798 + int OSUAnalysis::getPdgIdBinValue(int pdgId){
2799 +
2800 +  int binValue = -999;
2801 +  int absPdgId = fabs(pdgId);  
2802 +  if(pdgId == -1) binValue = 0;
2803 +  else if(absPdgId <= 6 ) binValue = absPdgId;
2804 +  else if(absPdgId == 11 ) binValue = 7;
2805 +  else if(absPdgId == 13 ) binValue = 8;
2806 +  else if(absPdgId == 15 ) binValue = 9;
2807 +  else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
2808 +  else if(absPdgId == 21 ) binValue = 11;
2809 +  else if(absPdgId == 22 ) binValue = 12;
2810 +  else if(absPdgId == 23 ) binValue = 13;
2811 +  else if(absPdgId == 24 ) binValue = 14;
2812 +  else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130  ) binValue = 15;
2813 +  else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400)  ) binValue = 16;
2814 +  else if(absPdgId > 400 && absPdgId < 500  ) binValue = 17;
2815 +  else if(absPdgId > 500 && absPdgId < 600  ) binValue = 18;
2816 +  else if(absPdgId > 1000 && absPdgId < 3000  ) binValue = 19;
2817 +  else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
2818 +  else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
2819 +  else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
2820 +
2821 +  else binValue = 23;
2822 +
2823 +  return binValue;
2824 +
2825 + }
2826  
2827  
2828   DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines