[ViewVC] Diff of: cvsroot/UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc

Comparing UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc (file contents):
Revision 1.12 by ahart, Fri Feb 8 11:07:51 2013 UTC vs.
Revision 1.32 by qpython, Tue Mar 19 18:05:23 2013 UTC

+  photons_ (cfg.getParameter<edm::InputTag> ("photons")),
+  superclusters_ (cfg.getParameter<edm::InputTag> ("superclusters")),
+  triggers_ (cfg.getParameter<edm::InputTag> ("triggers")),
-<
-<
->
+  puFile_ (cfg.getParameter<std::string> ("puFile")),
->
+  deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")),
->
+  muonSFFile_ (cfg.getParameter<std::string> ("muonSFFile")),
->
+  dataPU_ (cfg.getParameter<std::string> ("dataPU")),
->
+  electronSFID_ (cfg.getParameter<std::string> ("electronSFID")),
->
+  muonSF_ (cfg.getParameter<std::string> ("muonSF")),
->
+  dataset_ (cfg.getParameter<std::string> ("dataset")),
->
+  datasetType_ (cfg.getParameter<std::string> ("datasetType")),
+  channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels")),
-<
-<
+  histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets"))
-<
->
+  histogramSets_ (cfg.getParameter<vector<edm::ParameterSet> >("histogramSets")),
->
+  plotAllObjectsInPassingEvents_ (cfg.getParameter<bool> ("plotAllObjectsInPassingEvents")),
->
+  doPileupReweighting_ (cfg.getParameter<bool> ("doPileupReweighting"))
-–
+  TH1::SetDefaultSumw2 ();
-<
+  TH2::SetDefaultSumw2 ();
->
->
+  //create pile-up reweighting object, if necessary
->
+  if(datasetType_ != "data") {
->
+    if(doPileupReweighting_) puWeight_ = new PUWeight (puFile_, dataPU_, dataset_);
->
+    muonSFWeight_ = new MuonSFWeight (muonSFFile_, muonSF_);
->
+    electronSFWeight_ = new ElectronSFWeight ("53X", electronSFID_);
->
+  }
+  // Construct Cutflow Objects. These store the results of cut decisions and
+  // handle filling cut flow histograms.
+  std::vector<TFileDirectory> directories;
+  //always get vertex collection so we can assign the primary vertex in the event
-<
+  allNecessaryObjects.push_back("primaryvertexs");
->
+  objectsToGet.push_back("primaryvertexs");
->
->
+  //always make the plot of number of primary verticex (to check pile-up reweighting)
->
+  objectsToPlot.push_back("primaryvertexs");
-+
+  //always get the MC particles to do GEN-matching
-+
+  objectsToGet.push_back("mcparticles");
-+
-+
+  //always get the event collection to do pile-up reweighting
-+
+  objectsToGet.push_back("events");
+  //parse the histogram definitions
+  for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
-<
+    string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string>("inputCollection");
-<
+    objectsToPlot.push_back(tempInputCollection);
-<
+    allNecessaryObjects.push_back(tempInputCollection);
->
+    string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
->
+    if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
->
+    if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
->
+      objectsToGet.push_back(tempInputCollection);
->
+      objectsToPlot.push_back(tempInputCollection);
->
+      objectsToCut.push_back(tempInputCollection);
->
+    }
->
+    else{//pair of objects, need to add them both to the things to objectsToGet
->
+      int dashIndex = tempInputCollection.find("-");
->
+      int spaceIndex = tempInputCollection.find(" ");
->
+      int secondWordLength = spaceIndex - dashIndex;
->
+      objectsToGet.push_back(tempInputCollection);
->
+      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
->
+      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
->
+      objectsToPlot.push_back(tempInputCollection);
->
+      objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
->
+      objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
->
+      objectsToCut.push_back(tempInputCollection);
->
+      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
->
+      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
->
->
+      }
->
+    vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
-<
->
+    for(uint currentHistogram = 0; currentHistogram != histogramList_.size(); currentHistogram++){
+      histogram tempHistogram;
+      tempHistogram.title = histogramList_.at(currentHistogram).getParameter<string>("title");
+      tempHistogram.bins = histogramList_.at(currentHistogram).getParameter<vector<double> >("bins");
+      tempHistogram.inputVariables = histogramList_.at(currentHistogram).getParameter<vector<string> >("inputVariables");
-<
+      if(histogramList_.at(currentHistogram).exists("function"))
-<
+        tempHistogram.function = histogramList_.at(currentHistogram).getParameter<string>("function");
-<
+      else
-<
+        tempHistogram.function = "";
-<
->
+      histograms.push_back(tempHistogram);
-+
+  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
-+
+  sort( objectsToPlot.begin(), objectsToPlot.end() );
-+
+  objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
-+
-+
-+
-+
+  //add histograms with the gen-matched id, mother id, and grandmother id
-+
+  for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
-+
-+
+    string currentObject = objectsToPlot.at(currentObjectIndex);
-+
+    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
-+
-+
+    histogram tempIdHisto;
-+
+    histogram tempMomIdHisto;
-+
+    histogram tempGmaIdHisto;
-+
-+
+    tempIdHisto.inputCollection = currentObject;
-+
+    tempMomIdHisto.inputCollection = currentObject;
-+
+    tempGmaIdHisto.inputCollection = currentObject;
-+
-+
+    currentObject = currentObject.substr(0, currentObject.size()-1);
-+
+    tempIdHisto.name = currentObject+"GenMatchId";
-+
+    tempMomIdHisto.name = currentObject+"GenMatchMotherId";
-+
+    tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
-+
-+
+    currentObject.at(0) = toupper(currentObject.at(0));
-+
+    tempIdHisto.title = currentObject+" Gen-matched Particle";
-+
+    tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
-+
+    tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
-+
-+
+    int maxNum = 24;
-+
+    vector<double> binVector;
-+
+    binVector.push_back(maxNum);
-+
+    binVector.push_back(0);
-+
+    binVector.push_back(maxNum);
-+
-+
+    tempIdHisto.bins = binVector;
-+
+    tempIdHisto.inputVariables.push_back("genMatchedId");
-+
+    tempMomIdHisto.bins = binVector;
-+
+    tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
-+
+    tempGmaIdHisto.bins = binVector;
-+
+    tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
-+
-+
+    histograms.push_back(tempIdHisto);
-+
+    histograms.push_back(tempMomIdHisto);
-+
+    histograms.push_back(tempGmaIdHisto);
-+
-+
+  }
+      triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
+      for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
+        tempChannel.triggers.push_back(triggerNames.at(trigger));
-<
+      allNecessaryObjects.push_back("triggers");
->
+      objectsToGet.push_back("triggers");
-+
+    //create cutFlow for this channel
+    cutFlows_.push_back (new CutFlow (fs_, channelName));
+    std::map<std::string, TH2D*> twoDhistoMap;
+    twoDHists_.push_back(twoDhistoMap);
-+
+    //book all histograms included in the configuration
+    for(uint currentHistogramIndex = 0; currentHistogramIndex != histograms.size(); currentHistogramIndex++){
+      histogram currentHistogram = histograms.at(currentHistogramIndex);
+      int numBinsElements = currentHistogram.bins.size();
+      int numInputVariables = currentHistogram.inputVariables.size();
-+
+      if(numBinsElements != 3 && numBinsElements !=6) {
-<
+        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
-<
+        exit(0);
->
+        std::cout << "Error: Didn't find correct number of bin specifications for histogram named '" << currentHistogram.name << "'\n";
->
+        exit(0);
+      else if((numBinsElements == 3 && numInputVariables !=1) || (numBinsElements == 6 && numInputVariables!=2)){
-<
+        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
-<
+        exit(0);
->
+        std::cout << "Error: Didn't find correct number of input variables for histogram named '" << currentHistogram.name << "'\n";
->
+        exit(0);
-<
+      else if(numBinsElements == 3)
-<
+        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));
-<
+      else if(numBinsElements == 6)
-<
+        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));
->
+      else if(numBinsElements == 3){
->
+        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));
->
+      }
->
+      else if(numBinsElements == 6){
->
+        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));
->
->
+      }
->
->
->
+      if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
-+
+// bin      particle type
-+
+// ---      -------------
-+
+//  0        unmatched
-+
+//  1        u
-+
+//  2        d
-+
+//  3        s
-+
+//  4        c
-+
+//  5        b
-+
+//  6        t
-+
+//  7        e
-+
+//  8        mu
-+
+//  9        tau
-+
+// 10        nu
-+
+// 11        g
-+
+// 12        gamma
-+
+// 13        Z
-+
+// 14        W
-+
+// 15        light meson
-+
+// 16        K meson
-+
+// 17        D meson
-+
+// 18        B meson
-+
+// 19        light baryon
-+
+// 20        strange baryon
-+
+// 21        charm baryon
-+
+// 22        bottom baryon
-+
+// 23        other
-+
-+
+      vector<TString> labelArray;
-+
+      labelArray.push_back("unmatched");
-+
+      labelArray.push_back("u");
-+
+      labelArray.push_back("d");
-+
+      labelArray.push_back("s");
-+
+      labelArray.push_back("c");
-+
+      labelArray.push_back("b");
-+
+      labelArray.push_back("t");
-+
+      labelArray.push_back("e");
-+
+      labelArray.push_back("#mu");
-+
+      labelArray.push_back("#tau");
-+
+      labelArray.push_back("#nu");
-+
+      labelArray.push_back("g");
-+
+      labelArray.push_back("#gamma");
-+
+      labelArray.push_back("Z");
-+
+      labelArray.push_back("W");
-+
+      labelArray.push_back("light meson");
-+
+      labelArray.push_back("K meson");
-+
+      labelArray.push_back("D meson");
-+
+      labelArray.push_back("B meson");
-+
+      labelArray.push_back("light baryon");
-+
+      labelArray.push_back("strange baryon");
-+
+      labelArray.push_back("charm baryon");
-+
+      labelArray.push_back("bottom baryon");
-+
+      labelArray.push_back("other");
-+
+      for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
-+
+        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
-+
+      }
-+
+    //book a histogram for the number of each object type to be plotted
-+
+    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
+      string currentObject = objectsToPlot.at(currentObjectIndex);
+      int maxNum = 10;
+      if(currentObject == "mcparticles") maxNum = 50;
-<
+      currentObject[0] = toupper(currentObject[0]);
->
+      else if(currentObject == "primaryvertexs") maxNum = 50;
->
+      else if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
->
+      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
->
+      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
->
->
+      currentObject.at(0) = toupper(currentObject.at(0));
+      string histoName = "num" + currentObject;
-<
+      oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
->
->
+      if(histoName == "numPrimaryvertexs"){
->
+        string newHistoName = histoName + "BeforePileupCorrection";
->
+        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+" Before Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
->
+        newHistoName = histoName + "AfterPileupCorrection";
->
+        oneDHists_.at(currentChannel)[newHistoName] = directories.at(currentChannel).make<TH1D> (TString(newHistoName),channelLabel+" channel: Number of Selected "+currentObject+ " After Pileup Correction; # "+currentObject, maxNum, 0, maxNum);
->
+      }
->
+      else
->
+        oneDHists_.at(currentChannel)[histoName] = directories.at(currentChannel).make<TH1D> (TString(histoName),channelLabel+" channel: Number of Selected "+currentObject+"; # "+currentObject, maxNum, 0, maxNum);
-+
-+
-+
+    //get list of cuts for this channel
+    vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
+    //loop over and parse all cuts
+    for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
-–
+      cut tempCut;
-<
+     //store input collection for cut
-<
+      string inputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
-<
+      tempCut.inputCollection = inputCollection;
-<
+      allNecessaryObjects.push_back(inputCollection);
->
+      //store input collection for cut
->
+      string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
->
+      tempCut.inputCollection = tempInputCollection;
->
+      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
->
+        objectsToGet.push_back(tempInputCollection);
->
+        objectsToCut.push_back(tempInputCollection);
->
+      }
->
+      else{//pair of objects, need to add them both to the things to objectsToGet
->
+        int dashIndex = tempInputCollection.find("-");
->
+        int spaceIndex = tempInputCollection.find(" ");
->
+        int secondWordLength = spaceIndex - dashIndex;
->
+        objectsToGet.push_back(tempInputCollection);
->
+        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
->
+        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
->
+        objectsToCut.push_back(tempInputCollection);
->
+        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
->
+        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
->
->
+      }
->
->
+      //split cut string into parts and store them
+      string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
+      std::vector<string> cutStringVector = splitString(cutString);
-<
+      if(cutStringVector.size()!=3){
-<
+        std::cout << "Error: Didn't find three elements in the following cut string: '" <<cutString << "'\n";
-<
+        exit(0);
->
+      if(cutStringVector.size()!=3 && cutStringVector.size() % 4 !=3){
->
+        std::cout << "Error: Didn't find the expected number elements in the following cut string: '" << cutString << "'\n";
->
+        exit(0);
->
+      }
->
+      tempCut.numSubcuts = (cutStringVector.size()+1)/4;
->
+      for(int subcutIndex = 0; subcutIndex != tempCut.numSubcuts; subcutIndex++){//loop over all the pieces of the cut combined using &,|
->
+        int indexOffset = 4 * subcutIndex;
->
+        string currentVariableString = cutStringVector.at(indexOffset);
->
+        if(currentVariableString.find("(")==std::string::npos){
->
+          tempCut.functions.push_back("");//no function was specified
->
+          tempCut.variables.push_back(currentVariableString);// variable to cut on
->
+        }
->
+        else{
->
+          tempCut.functions.push_back(currentVariableString.substr(0,currentVariableString.find("(")));//function comes before the "("
->
+          string tempVariable = currentVariableString.substr(currentVariableString.find("(")+1);//get rest of string
->
+          tempCut.variables.push_back(tempVariable.substr(0,tempVariable.size()-1));//remove trailing ")"
->
+        }
->
+        tempCut.comparativeOperators.push_back(cutStringVector.at(indexOffset+1));// comparison to make
->
+        tempCut.cutValues.push_back(atof(cutStringVector.at(indexOffset+2).c_str()));// threshold value to pass cut
->
+        if(subcutIndex != 0) tempCut.logicalOperators.push_back(cutStringVector.at(indexOffset-1)); // logical comparison (and, or)
-–
-–
+      tempCut.variable = cutStringVector.at(0);// variable to cut on
-–
+      tempCut.comparativeOperator = cutStringVector.at(1);// comparison to make
-–
+      tempCut.cutValue = atof(cutStringVector.at(2).c_str());// threshold value to pass cut
+      //get number of objects required to pass cut for event to pass
+      string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
+      std::vector<string> numberRequiredVector = splitString(numberRequiredString);
+      if(numberRequiredVector.size()!=2){
-<
+        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
->
+        std::cout << "Error: Didn't find two elements in the following number requirement string: '" << numberRequiredString << "'\n";
+        exit(0);
-–
+      // determine number required if comparison contains "="
+      int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
-–
+      if(numberRequiredVector.at(0) == ">") numberRequiredInt++;
-–
+      else if(numberRequiredVector.at(0) == "<") numberRequiredInt--;
-–
+      tempCut.numberRequired = numberRequiredInt;// number of objects required to pass the cut
+      tempCut.eventComparativeOperator = numberRequiredVector.at(0);// comparison to make
-<
+      if(cuts_.at(currentCut).exists("function")){
-<
+        tempCut.function = cuts_.at(currentCut).getParameter<string> ("function");
-<
+      }
-<
+      else tempCut.function = "";
->
+      string tempCutName;
+      if(cuts_.at(currentCut).exists("alias")){
+        tempCutName = cuts_.at(currentCut).getParameter<string> ("alias");
+      else{
+        //construct string for cutflow table
+        bool plural = numberRequiredInt != 1;
-<
+        string collectionString = plural ? inputCollection : inputCollection.substr(0, inputCollection.size()-1);
->
+        string collectionString = plural ? tempInputCollection : tempInputCollection.substr(0, tempInputCollection.size()-1);
+        string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
+        tempCutName = cutName;
+    }//end loop over cuts
-–
+    channels.push_back(tempChannel);
+    tempChannel.cuts.clear();
+  }//end loop over channels
-<
-<
-<
+  //make unique vector of objects we need to cut on (so we make sure to get them from the event)
-<
+  sort( allNecessaryObjects.begin(), allNecessaryObjects.end() );
-<
+  allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() );
->
+  //make unique vector of objects we need to get from the event
->
+  sort( objectsToGet.begin(), objectsToGet.end() );
->
+  objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
->
+  //make unique vector of objects we need to cut on
->
+  sort( objectsToCut.begin(), objectsToCut.end() );
->
+  objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
+OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
-–
+  // Retrieve necessary collections from the event.
-<
+  edm::Handle<BNtriggerCollection> triggers;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "triggers") != allNecessaryObjects.end())
->
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
+    event.getByLabel (triggers_, triggers);
-<
+  edm::Handle<BNjetCollection> jets;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
+    event.getByLabel (jets_, jets);
-<
+  edm::Handle<BNmuonCollection> muons;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
+    event.getByLabel (muons_, muons);
-<
+  edm::Handle<BNelectronCollection> electrons;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
+    event.getByLabel (electrons_, electrons);
-<
+  edm::Handle<BNeventCollection> events;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
+    event.getByLabel (events_, events);
-<
+  edm::Handle<BNtauCollection> taus;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
+    event.getByLabel (taus_, taus);
-<
+  edm::Handle<BNmetCollection> mets;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
+    event.getByLabel (mets_, mets);
-<
+  edm::Handle<BNtrackCollection> tracks;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
+    event.getByLabel (tracks_, tracks);
-<
+  edm::Handle<BNgenjetCollection> genjets;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
+    event.getByLabel (genjets_, genjets);
-<
+  edm::Handle<BNmcparticleCollection> mcparticles;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
+    event.getByLabel (mcparticles_, mcparticles);
-<
+  edm::Handle<BNprimaryvertexCollection> primaryvertexs;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
+    event.getByLabel (primaryvertexs_, primaryvertexs);
-<
+  edm::Handle<BNbxlumiCollection> bxlumis;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
+    event.getByLabel (bxlumis_, bxlumis);
-<
+  edm::Handle<BNphotonCollection> photons;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
+    event.getByLabel (photons_, photons);
-<
+  edm::Handle<BNsuperclusterCollection> superclusters;
-<
+  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end())
->
+  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
+    event.getByLabel (superclusters_, superclusters);
-+
+  //get pile-up event weight
-+
+  double scaleFactor = 1.00;
-+
+  if(doPileupReweighting_ && datasetType_ != "data")
-+
+    scaleFactor = puWeight_->at (events->at (0).numTruePV);
-+
+  //loop over all channels
+    channel currentChannel = channels.at(currentChannelIndex);
+    flagMap individualFlags;
-–
+    flagMap cumulativeFlags;
+    counterMap passingCounter;
-+
+    cumulativeFlags.clear ();
+    bool triggerDecision = true;
+    if(currentChannel.triggers.size() != 0){  //triggers specified
+    //loop over all cuts
-+
-+
-+
+    for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
+      cut currentCut = currentChannel.cuts.at(currentCutIndex);
-<
+      for(uint currentObjectIndex = 0; currentObjectIndex != allNecessaryObjects.size(); currentObjectIndex++){
->
+      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
->
+        string currentObject = objectsToCut.at(currentObjectIndex);
-<
+        string currentObject = allNecessaryObjects.at(currentObjectIndex);
->
+        //initialize maps to get ready to set cuts
+        individualFlags[currentObject].push_back (vector<bool> ());
+        cumulativeFlags[currentObject].push_back (vector<bool> ());
-+
+      }
-+
+      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
-+
+        string currentObject = objectsToCut.at(currentObjectIndex);
-+
-+
+        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
-+
-+
+        if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
+        else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
+        else if(currentObject == "electrons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),"electrons");
+        else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
-<
+      }
->
+        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
->
+                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
->
+                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
->
+                                                                   "muon-muon pairs");
->
+        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
->
+                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
->
+                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
->
+                                                                           "electron-electron pairs");
->
+        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
->
+                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
->
+                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
->
+                                                                       "electron-muon pairs");
-+
+      }
-–
+    }//end loop over all cuts
-+
+    }//end loop over all cuts
+    //use cumulative flags to apply cuts at event level
+    eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
-+
+    for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
+      //loop over all objects and count how many passed the cumulative selection up to this point
+      cut currentCut = currentChannel.cuts.at(currentCutIndex);
+      int numberPassing = 0;
-<
+      for (uint object = 0; object != cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).size() ; object++)
-<
+          if(cumulativeFlags[currentCut.inputCollection].at(currentCutIndex).at(object)) numberPassing++;
-<
->
+      for (uint object = 0; object != cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).size() ; object++){
->
+          if(cumulativeFlags.at(currentCut.inputCollection).at(currentCutIndex).at(object)) numberPassing++;
->
+      }
+      bool cutDecision = evaluateComparison(numberPassing,currentCut.eventComparativeOperator,currentCut.numberRequired);
+      cutFlows_.at(currentChannelIndex)->at (currentCut.name) = cutDecision;
-<
+    cutFlows_.at(currentChannelIndex)->fillCutFlow();
-<
->
+    cutFlows_.at(currentChannelIndex)->fillCutFlow(scaleFactor);
-–
+    if(!eventPassedAllCuts)continue;
-<
+    //set position of primary vertex in event, in order to calculate quantities relative to it
-<
+    primaryVertex_ = 0;
-<
+    vector<bool> vertexFlags = cumulativeFlags["primaryvertexs"].back();
-<
+    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
-<
+      if(!vertexFlags.at(vertexIndex)) continue;
-<
+      primaryVertex_ = new BNprimaryvertex (primaryvertexs->at (vertexIndex));
-<
+      break;
-<
+    }
->
+    if(!eventPassedAllCuts)continue;
-+
+    //if(datasetType_ != "data") {
-+
+    //  scaleFactor *= muonSFWeight_->at (chosenMuon ()->eta);
-+
+    //  scaleFactor *= electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
-+
+    //}
+    //filling histograms
+    for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
+      histogram currentHistogram = histograms.at(histogramIndex);
+      if(currentHistogram.inputVariables.size() == 1){
-<
+        TH1D* histo;
-<
+        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
-<
+        if(currentHistogram.inputCollection == "jets") fillHistogram(histo,currentHistogram,jets.product());
-<
+        else if(currentHistogram.inputCollection == "muons") fillHistogram(histo,currentHistogram,muons.product());
-<
+        else if(currentHistogram.inputCollection == "electrons") fillHistogram(histo,currentHistogram,electrons.product());
-<
+        else if(currentHistogram.inputCollection == "events") fillHistogram(histo,currentHistogram,events.product());
-<
+        else if(currentHistogram.inputCollection == "taus") fillHistogram(histo,currentHistogram,taus.product());
-<
+        else if(currentHistogram.inputCollection == "mets") fillHistogram(histo,currentHistogram,mets.product());
-<
+        else if(currentHistogram.inputCollection == "tracks") fillHistogram(histo,currentHistogram,tracks.product());
-<
+        else if(currentHistogram.inputCollection == "genjets") fillHistogram(histo,currentHistogram,genjets.product());
-<
+        else if(currentHistogram.inputCollection == "mcparticles") fillHistogram(histo,currentHistogram,mcparticles.product());
-<
+        else if(currentHistogram.inputCollection == "primaryvertexs") fillHistogram(histo,currentHistogram,primaryvertexs.product());
-<
+        else if(currentHistogram.inputCollection == "bxlumis") fillHistogram(histo,currentHistogram,bxlumis.product());
-<
+        else if(currentHistogram.inputCollection == "photons") fillHistogram(histo,currentHistogram,photons.product());
-<
+        else if(currentHistogram.inputCollection == "superclusters") fillHistogram(histo,currentHistogram,superclusters.product());
->
+        TH1D* histo;
->
+        histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
->
->
->
->
+        if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
->
+         else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
->
+                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
->
+                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\
->
+                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\
->
+                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "electron-muon pairs") fill1DHistogram(histo,currentHistogram, electrons.product(),muons.product(), \
->
+                                                                                              cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(),
->
+                                                                                              cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "events") fill1DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "taus") fill1DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "mets") fill1DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "tracks") fill1DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "genjets") fill1DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "mcparticles") fill1DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "primaryvertexs") fill1DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "bxlumis") fill1DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "photons") fill1DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
->
+        else if(currentHistogram.inputCollection == "superclusters") fill1DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
+      else if(currentHistogram.inputVariables.size() == 2){
-<
+        TH2D* histo;
-<
+        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
-<
+        if(currentHistogram.inputCollection == "jets") fillHistogram(histo,currentHistogram,jets.product());
-<
+        else if(currentHistogram.inputCollection == "muons") fillHistogram(histo,currentHistogram,muons.product());
-<
+        else if(currentHistogram.inputCollection == "electrons") fillHistogram(histo,currentHistogram,electrons.product());
-<
+        else if(currentHistogram.inputCollection == "events") fillHistogram(histo,currentHistogram,events.product());
-<
+        else if(currentHistogram.inputCollection == "taus") fillHistogram(histo,currentHistogram,taus.product());
-<
+        else if(currentHistogram.inputCollection == "mets") fillHistogram(histo,currentHistogram,mets.product());
-<
+        else if(currentHistogram.inputCollection == "tracks") fillHistogram(histo,currentHistogram,tracks.product());
-<
+        else if(currentHistogram.inputCollection == "genjets") fillHistogram(histo,currentHistogram,genjets.product());
-<
+        else if(currentHistogram.inputCollection == "mcparticles") fillHistogram(histo,currentHistogram,mcparticles.product());
-<
+        else if(currentHistogram.inputCollection == "primaryvertexs") fillHistogram(histo,currentHistogram,primaryvertexs.product());
-<
+        else if(currentHistogram.inputCollection == "bxlumis") fillHistogram(histo,currentHistogram,bxlumis.product());
-<
+        else if(currentHistogram.inputCollection == "photons") fillHistogram(histo,currentHistogram,photons.product());
-<
+        else if(currentHistogram.inputCollection == "superclusters") fillHistogram(histo,currentHistogram,superclusters.product());
-<
+      }
-<
+    }
->
+        TH2D* histo;
->
+        histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
-+
-+
+        if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
-+
+                                                                                       cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \
-+
+                                                                                       cumulativeFlags.at("muon-muon pairs").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "electrons") fill2DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "electron-electron pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),electrons.product(), \
-+
+                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(), \
-+
+                                                                                               cumulativeFlags.at("electron-electron pairs").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "electron-muon pairs") fill2DHistogram(histo,currentHistogram,electrons.product(),muons.product(), \
-+
+                                                                                               cumulativeFlags.at("electrons").back(),cumulativeFlags.at("muons").back(), \
-+
+                                                                                               cumulativeFlags.at("electron-muon pairs").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "events") fill2DHistogram(histo,currentHistogram,events.product(),cumulativeFlags.at("events").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "taus") fill2DHistogram(histo,currentHistogram,taus.product(),cumulativeFlags.at("taus").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "mets") fill2DHistogram(histo,currentHistogram,mets.product(),cumulativeFlags.at("mets").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "tracks") fill2DHistogram(histo,currentHistogram,tracks.product(),cumulativeFlags.at("tracks").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "genjets") fill2DHistogram(histo,currentHistogram,genjets.product(),cumulativeFlags.at("genjets").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "mcparticles") fill2DHistogram(histo,currentHistogram,mcparticles.product(),cumulativeFlags.at("mcparticles").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "primaryvertexs") fill2DHistogram(histo,currentHistogram,primaryvertexs.product(),cumulativeFlags.at("primaryvertexs").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "bxlumis") fill2DHistogram(histo,currentHistogram,bxlumis.product(),cumulativeFlags.at("bxlumis").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "photons") fill2DHistogram(histo,currentHistogram,photons.product(),cumulativeFlags.at("photons").back(),scaleFactor);
-+
+        else if(currentHistogram.inputCollection == "superclusters") fill2DHistogram(histo,currentHistogram,superclusters.product(),cumulativeFlags.at("superclusters").back(),scaleFactor);
-+
+      }
-+
+    }
-+
+    //fills histograms with the sizes of collections
+    for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
+      string currentObject = objectsToPlot.at(currentObjectIndex);
-<
+      string tempCurrentObject = currentObject;
-<
+      tempCurrentObject[0] = toupper(tempCurrentObject[0]);
->
->
+      string objectToPlot = "";
->
->
+      if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
->
+      else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
->
+      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
->
+      else objectToPlot = currentObject;
->
+      string tempCurrentObject = objectToPlot;
->
+      tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
+      string histoName = "num" + tempCurrentObject;
-<
+      if(currentObject == "jets") oneDHists_.at(currentChannelIndex)[histoName]->Fill(jets->size());
-<
+      else if(currentObject == "muons") oneDHists_.at(currentChannelIndex)[histoName]->Fill(muons->size());
-<
+      else if(currentObject == "electrons") oneDHists_.at(currentChannelIndex)[histoName]->Fill(electrons->size());
-<
+      else if(currentObject == "events") oneDHists_.at(currentChannelIndex)[histoName]->Fill(events->size());
-<
+      else if(currentObject == "taus") oneDHists_.at(currentChannelIndex)[histoName]->Fill(taus->size());
-<
+      else if(currentObject == "mets") oneDHists_.at(currentChannelIndex)[histoName]->Fill(mets->size());
-<
+      else if(currentObject == "tracks") oneDHists_.at(currentChannelIndex)[histoName]->Fill(tracks->size());
-<
+      else if(currentObject == "genjets") oneDHists_.at(currentChannelIndex)[histoName]->Fill(genjets->size());
-<
+      else if(currentObject == "mcparticles") oneDHists_.at(currentChannelIndex)[histoName]->Fill(mcparticles->size());
-<
+      else if(currentObject == "primaryvertexs") oneDHists_.at(currentChannelIndex)[histoName]->Fill(primaryvertexs->size());
-<
+      else if(currentObject == "bxlumis") oneDHists_.at(currentChannelIndex)[histoName]->Fill(bxlumis->size());
-<
+      else if(currentObject == "photons") oneDHists_.at(currentChannelIndex)[histoName]->Fill(photons->size());
-<
+      else if(currentObject == "superclusters") oneDHists_.at(currentChannelIndex)[histoName]->Fill(superclusters->size());
->
->
->
->
+      //set position of primary vertex in event, in order to calculate quantities relative to it
->
+      if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
->
+        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
->
+        int numToPlot = 0;
->
+        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
->
+          if(lastCutFlags.at(currentFlag)) numToPlot++;
->
+        }
->
+        if(objectToPlot == "primaryvertexs"){
->
+          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
->
+          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
->
+        }
->
+        else
->
+          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,scaleFactor);
->
+      }
->
+      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),scaleFactor);
->
+      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),scaleFactor);
->
+      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,scaleFactor);
->
+      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),scaleFactor);
->
+      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,scaleFactor);
->
+      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),scaleFactor);
->
+      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),scaleFactor);
->
+      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),scaleFactor);
->
+      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),scaleFactor);
->
+      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),scaleFactor);
->
+      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),scaleFactor);
->
+      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),scaleFactor);
->
+      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),scaleFactor);
->
+      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),scaleFactor);
->
+      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),scaleFactor);
->
+      else if(objectToPlot == "primaryvertexs"){
->
+        oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
->
+        oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),scaleFactor);
->
+      }
-+
-+
-+
-+
+  } //end loop over channel
-<
+  masterCutFlow_->fillCutFlow();
->
+  masterCutFlow_->fillCutFlow(scaleFactor);
->
+  else if(comparison == "<")  return testValue <  cutValue;
+  else if(comparison == "<=") return testValue <= cutValue;
+  else if(comparison == "==") return testValue == cutValue;
-+
+  else if(comparison == "=") return testValue == cutValue;
+  else if(comparison == "!=") return testValue != cutValue;
+  else {std::cout << "WARNING: invalid comparison operator '" << comparison << "'\n"; return false;}
+  else if(variable == "ecalIso") value = object->ecalIso;
+  else if(variable == "hcalIso") value = object->hcalIso;
+  else if(variable == "caloIso") value = object->caloIso;
-<
+  else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
->
+  else if(variable == "trackIsDR03") value = object->trackIsoDR03;
+  else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
+  else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
+  else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
+  else if(variable == "time_ndof") value = object->time_ndof;
+  //user-defined variables
-<
+  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
-<
+  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
-<
+  else if(variable == "detIso") value = (object->trackIso + object->caloIso) / object->pt;
-<
+  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
->
+  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
->
+  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
->
+  else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
->
+  else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
+  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
-+
+  else if(variable == "metMT") {
-+
+    const BNmet *met = chosenMET ();
-+
+    if (met)
-+
+      {
-+
+        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
-+
+                       p2 (met->px, met->py, 0.0, met->pt);
-+
-+
+        value = (p1 + p2).Mt ();
-+
+      }
-+
+    else
-+
+      value = -999;
-+
+  }
-+
-+
-+
-+
+  else if(variable == "correctedD0VertexInEBPlus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBPlus"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEPlus"){
-+
+    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBPlus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBPlus"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEPlus"){
-+
+    if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBMinus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBMinus"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEMinus"){
-+
+    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBMinus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBMinus"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEMinus"){
-+
+    if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
-+
+  else if(variable == "correctedD0VertexInEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEPositiveCharge"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEENegativeCharge"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEENegativeCharge"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "tightID") {
-<
+    value = object->isGlobalMuon > 0            \
-<
+      && object->isPFMuon > 0                   \
-<
+      && object->normalizedChi2 < 10            \
-<
+      && object->numberOfValidMuonHits > 0      \
-<
+      && object->numberOfMatchedStations > 1    \
-<
+      && fabs(object->correctedD0Vertex) < 0.2  \
-<
+      && fabs(object->correctedDZ) < 0.5        \
-<
+      && object->numberOfValidPixelHits > 0     \
->
+    value = object->isGlobalMuon > 0                \
->
+      && object->isPFMuon > 0                        \
->
+      && object->normalizedChi2 < 10                \
->
+      && object->numberOfValidMuonHits > 0        \
->
+      && object->numberOfMatchedStations > 1        \
->
+      && fabs(object->correctedD0Vertex) < 0.2        \
->
+      && fabs(object->correctedDZ) < 0.5        \
->
+      && object->numberOfValidPixelHits > 0                \
+      && object->numberOfLayersWithMeasurement > 5;
-–
-–
+  else if(variable == "tightIDdisplaced"){
-<
+    value = object->isGlobalMuon > 0            \
-<
+      && object->isPFMuon > 0                   \
-<
+      && object->normalizedChi2 < 10            \
-<
+      && object->numberOfValidMuonHits > 0      \
-<
+      && object->numberOfMatchedStations > 1    \
-<
+      && object->numberOfValidPixelHits > 0     \
->
+    value = object->isGlobalMuon > 0                \
->
+      && object->normalizedChi2 < 10                \
->
+      && object->numberOfValidMuonHits > 0        \
->
+      && object->numberOfMatchedStations > 1        \
->
+      && object->numberOfValidPixelHits > 0        \
+      && object->numberOfLayersWithMeasurement > 5;
-+
-+
-+
+  else if(variable == "genMatchedId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).id);
-+
+  }
-+
+  else if(variable == "genMatchedMotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
-+
+  }
-+
+  else if(variable == "genMatchedGrandmotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else if(fabs(mcparticles->at(index).motherId) == 15){
-+
+      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
-+
+      if(motherIndex == -1) value = 0;
-+
+      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
-+
+    }
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
-+
+  }
-+
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  value = applyFunction(function, value);
+  else if(variable == "passConvVeto") value = object->passConvVeto;
+  //user-defined variables
-<
+  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
-<
+  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
-<
+  else if(variable == "detIso") value = (object->trackIso + object->caloIso) / object->pt;
->
+  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
->
+  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenVertex ()->xError, chosenVertex ()->yError));
->
+  else if(variable == "detIso") value = (object->trackIso) / object->pt;
+  else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
-+
+  else if(variable == "metMT") {
-+
+    const BNmet *met = chosenMET ();
-+
+    if (met)
-+
+      {
-+
+        TLorentzVector p1 (object->px, object->py, object->pz, object->energy),
-+
+                       p2 (met->px, met->py, 0.0, met->pt);
-+
-+
+        value = (p1 + p2).Mt ();
-+
+      }
-+
+    else
-+
+      value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){
-+
+    if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBPlus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBPlus"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEPlus"){
-+
+    if(object->isEE && object->eta > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBPlus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBPlus"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEPlus"){
-+
+    if(object->isEE && object->eta > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBMinus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBMinus"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEMinus"){
-+
+    if(object->isEE && object->eta < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBMinus"){
-+
+    if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBMinus"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEMinus"){
-+
+    if(object->isEE && object->eta < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "tightID"){
-+
+    if (object->isEB)
-+
+      {
-+
+        value = fabs(object->delEtaIn) < 0.004 \
-+
+          && fabs (object->delPhiIn) < 0.03 \
-+
+          && object->scSigmaIEtaIEta < 0.01 \
-+
+          && object->hadOverEm < 0.12 \
-+
+          && fabs (object->correctedD0Vertex) < 0.02 \
-+
+          && fabs (object->correctedDZ) < 0.1 \
-+
+          && object->absInvEMinusInvPin < 0.05 \
-+
+          && object->passConvVeto;
-+
+      }
-+
+    else
-+
+      {
-+
+        value = fabs(object->delEtaIn) < 0.005 \
-+
+          && fabs (object->delPhiIn) < 0.02 \
-+
+          && object->scSigmaIEtaIEta < 0.03 \
-+
+          && object->hadOverEm < 0.10 \
-+
+          && fabs (object->correctedD0Vertex) < 0.02 \
-+
+          && fabs (object->correctedDZ) < 0.1 \
-+
+          && object->absInvEMinusInvPin < 0.05 \
-+
+          && object->passConvVeto;
-+
+      }
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBPositiveCharge"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEEPositiveCharge"){
-+
+    if(object->isEE && object->charge > 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBPositiveCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEEPositiveCharge"){
-+
+    if(object->isEE && object->charge > 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0VertexInEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexOutEBNegativeCharge"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0VertexEENegativeCharge"){
-+
+    if(object->isEE && object->charge < 0) value = object->correctedD0Vertex;
-+
+    else value = -999;
-+
+  }
-+
-+
+  else if(variable == "correctedD0BeamspotInEBNegativeCharge"){
-+
+    if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){
-+
+    if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "correctedD0BeamspotEENegativeCharge"){
-+
+    if(object->isEE && object->charge < 0) value = object->correctedD0;
-+
+    else value = -999;
-+
+  }
-+
-+
-+
+  else if(variable == "tightIDdisplaced"){
-+
+    if (object->isEB)
-+
+      {
-+
+        value = fabs(object->delEtaIn) < 0.004 \
-+
+          && fabs (object->delPhiIn) < 0.03 \
-+
+          && object->scSigmaIEtaIEta < 0.01 \
-+
+          && object->hadOverEm < 0.12 \
-+
+          && object->absInvEMinusInvPin < 0.05;
-+
+      }
-+
+    else
-+
+      {
-+
+        value = fabs (object->delEtaIn) < 0.005 \
-+
+          && fabs (object->delPhiIn) < 0.02 \
-+
+          && object->scSigmaIEtaIEta < 0.03 \
-+
+          && object->hadOverEm < 0.10 \
-+
+          && object->absInvEMinusInvPin < 0.05;
-+
+      }
-+
+  }
-+
-+
+  else if(variable == "genMatchedId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).id);
-+
+  }
-+
+  else if(variable == "genMatchedMotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
-+
+  }
-+
+  else if(variable == "genMatchedGrandmotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else if(fabs(mcparticles->at(index).motherId) == 15){
-+
+      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
-+
+      if(motherIndex == -1) value = 0;
-+
+      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
-+
+    }
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
-+
+  }
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  else if(variable == "id1") value = object->id1;
+  else if(variable == "id2") value = object->id2;
+  else if(variable == "evt") value = object->evt;
-+
+  else if(variable == "puScaleFactor"){
-+
+    if(datasetType_ != "data")
-+
+      value = puWeight_->at (events->at (0).numTruePV);
-+
+    else
-+
+      value = 1.0;
-+
+  }
-+
+  else if(variable == "muonScaleFactor"){
-+
+    if(datasetType_ != "data")
-+
+      value = muonSFWeight_->at (chosenMuon ()->eta);
-+
+    else
-+
+      value = 1.0;
-+
+  }
-+
+  else if(variable == "electronScaleFactor"){
-+
+    if(datasetType_ != "data")
-+
+      value = electronSFWeight_->at (chosenElectron ()->eta, chosenElectron ()->pt);
-+
+    else
-+
+      value = 1.0;
-+
+  }
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
-+
+  else if(variable == "genMatchedId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).id);
-+
+  }
-+
+  else if(variable == "genMatchedMotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
-+
+  }
-+
+  else if(variable == "genMatchedGrandmotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else if(fabs(mcparticles->at(index).motherId) == 15){
-+
+      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
-+
+      if(motherIndex == -1) value = 0;
-+
+      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
-+
+    }
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
-+
+  }
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  value = applyFunction(function, value);
+OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){
+  double value = 0.0;
-+
+  double pMag = sqrt(object->pt * object->pt +
-+
+                         object->pz * object->pz);
+  if(variable == "pt") value = object->pt;
+  else if(variable == "px") value = object->px;
+  else if(variable == "isHighPurity") value = object->isHighPurity;
-+
+  //additional BNs info for disappTrks
-+
+  else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution;
-+
+  else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3;
-+
+  else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3;
-+
+  else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4;
-+
+  else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4;
-+
+  else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5;
-+
+  else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5;
-+
+  else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter;
-+
+  else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner;
-+
+  else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle;
-+
-+
-+
+  //user defined variables
-+
+  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;
-+
+  else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz;
-+
+  else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5);
-+
+  else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag);
-+
-+
-+
-+
-+
-+
+  else if(variable == "isIso") value = getTrkIsIso(object, tracks.product());
-+
+  else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object);
-+
+  else if(variable == "ptErrorByPt") value = (object->ptError/object->pt);
-+
+  else if(variable == "ptError") value = object->ptError;
-+
+  else if(variable == "ptRes") value = getTrkPtRes(object);
-+
-+
-+
+  else if(variable == "genMatchedId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).id);
-+
+  }
-+
+  else if(variable == "genMatchedMotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
-+
+  }
-+
+  else if(variable == "genMatchedGrandmotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else if(fabs(mcparticles->at(index).motherId) == 15){
-+
+      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
-+
+      if(motherIndex == -1) value = 0;
-+
+      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
-+
+    }
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
-+
+  }
-+
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  value = applyFunction(function, value);
+  //user-defined variables
+  else if (variable == "d0"){
-<
+    double vx = object->vx - primaryVertex_->x,
-<
+      vy = object->vy - primaryVertex_->y,
->
+    double vx = object->vx - chosenVertex ()->x,
->
+      vy = object->vy - chosenVertex ()->y,
+      px = object->px,
+      py = object->py,
+      pt = object->pt;
+    value = (-vx * py + vy * px) / pt;
+  else if (variable == "dz"){
-<
+    double vx = object->vx - primaryVertex_->x,
-<
+      vy = object->vy - primaryVertex_->y,
-<
+      vz = object->vz - primaryVertex_->z,
->
+    double vx = object->vx - chosenVertex ()->x,
->
+      vy = object->vy - chosenVertex ()->y,
->
+      vz = object->vz - chosenVertex ()->z,
+      px = object->px,
+      py = object->py,
+      pz = object->pz,
+      pt = object->pt;
+    value = vz - (vx * px + vy * py)/pt * (pz/pt);
-<
->
+  else if(variable == "v0"){
->
+    value = sqrt(object->vx*object->vx + object->vy*object->vy);
->
+  }
->
+  else if(variable == "deltaV0"){
->
+    value = sqrt((object->vx-chosenVertex ()->x)*(object->vx-chosenVertex ()->x) + (object->vy-chosenVertex ()->y)*(object->vy-chosenVertex ()->y));
->
+  }
->
+  else if (variable == "deltaVx"){
->
+    value = object->vx - chosenVertex ()->x;
->
+  }
->
+  else if (variable == "deltaVy"){
->
+    value = object->vy - chosenVertex ()->y;
->
+  }
->
+  else if (variable == "deltaVz"){
->
+    value = object->vz - chosenVertex ()->z;
->
+  }
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
-+
+  else if(variable == "genMatchedId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).id);
-+
+  }
-+
+  else if(variable == "genMatchedMotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
-+
+  }
-+
+  else if(variable == "genMatchedGrandmotherId"){
-+
+    int index = getGenMatchedParticleIndex(object);
-+
+    if(index == -1) value = 0;
-+
+    else if(fabs(mcparticles->at(index).motherId) == 15){
-+
+      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
-+
+      if(motherIndex == -1) value = 0;
-+
+      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
-+
+    }
-+
+    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
-+
+  }
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
+  value = applyFunction(function, value);
+double
-+
+OSUAnalysis::valueLookup (const BNmuon* object1, const BNmuon* object2, string variable, string function){
-+
-+
+  double value = 0.0;
-+
-+
+  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
-+
+  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
-+
+  else if(variable == "invMass"){
-+
+    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
-+
+    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
-+
+    value = (fourVector1 + fourVector2).M();
-+
+  }
-+
+  else if(variable == "threeDAngle")
-+
+    {
-+
+      TVector3 threeVector1(object1->px, object1->py, object1->pz);
-+
+      TVector3 threeVector2(object2->px, object2->py, object2->pz);
-+
+      value = (threeVector1.Angle(threeVector2));
-+
+    }
-+
+  else if(variable == "alpha")
-+
+    {
-+
+      static const double pi = 3.1415926535897932384626433832795028841971693993751058;
-+
+      TVector3 threeVector1(object1->px, object1->py, object1->pz);
-+
+      TVector3 threeVector2(object2->px, object2->py, object2->pz);
-+
+      value = (pi-threeVector1.Angle(threeVector2));
-+
+    }
-+
+  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
-+
+  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
-+
+  else if(variable == "d0Sign"){
-+
+    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
-+
+    if(d0Sign < 0) value = -0.5;
-+
+    else if (d0Sign > 0) value = 0.5;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "chargeProduct"){
-+
+    value = object1->charge*object2->charge;
-+
+  }
-+
+  else if(variable == "muon1CorrectedD0Vertex"){
-+
+    value = object1->correctedD0Vertex;
-+
+  }
-+
+  else if(variable == "muon2CorrectedD0Vertex"){
-+
+    value = object2->correctedD0Vertex;
-+
+  }
-+
+ else if(variable == "muon1timeAtIpInOut"){
-+
+    value = object1->timeAtIpInOut;
-+
+  }
-+
+ else if(variable == "muon2timeAtIpInOut"){
-+
+    value = object2->timeAtIpInOut;
-+
+  }
-+
+ else if(variable == "muon1correctedD0")
-+
+   {
-+
+     value = object1->correctedD0;
-+
+   }
-+
+ else if(variable == "muon1correctedD0Vertex")
-+
+   {
-+
+     value = object1->correctedD0Vertex;
-+
+   }
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
-+
-+
+  value = applyFunction(function, value);
-+
-+
+  return value;
-+
+}
-+
-+
+double
-+
+OSUAnalysis::valueLookup (const BNelectron* object1, const BNelectron* object2, string variable, string function){
-+
-+
+  double value = 0.0;
-+
-+
+  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
-+
+  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
-+
+  else if(variable == "invMass"){
-+
+    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
-+
+    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
-+
+    value = (fourVector1 + fourVector2).M();
-+
+  }
-+
+  else if(variable == "threeDAngle")
-+
+    {
-+
+      TVector3 threeVector1(object1->px, object1->py, object1->pz);
-+
+      TVector3 threeVector2(object2->px, object2->py, object2->pz);
-+
+      value = (threeVector1.Angle(threeVector2));
-+
+    }
-+
+  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
-+
+  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
-+
+  else if(variable == "d0Sign"){
-+
+    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
-+
+    if(d0Sign < 0) value = -0.5;
-+
+    else if (d0Sign > 0) value = 0.5;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "chargeProduct"){
-+
+    value = object1->charge*object2->charge;
-+
+  }
-+
+  else if(variable == "electron1CorrectedD0Vertex"){
-+
+    value = object1->correctedD0Vertex;
-+
+  }
-+
+  else if(variable == "electron2CorrectedD0Vertex"){
-+
+    value = object2->correctedD0Vertex;
-+
+  }
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
-+
-+
+  value = applyFunction(function, value);
-+
-+
+  return value;
-+
+}
-+
-+
+double
-+
+OSUAnalysis::valueLookup (const BNelectron* object1, const BNmuon* object2, string variable, string function){
-+
-+
+  double value = 0.0;
-+
-+
+  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
-+
+  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
-+
+  else if(variable == "invMass"){
-+
+    TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
-+
+    TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
-+
+    value = (fourVector1 + fourVector2).M();
-+
+  }
-+
+  else if(variable == "threeDAngle")
-+
+    {
-+
+      TVector3 threeVector1(object1->px, object1->py, object1->pz);
-+
+      TVector3 threeVector2(object2->px, object2->py, object2->pz);
-+
+      value = (threeVector1.Angle(threeVector2));
-+
+    }
-+
+  else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
-+
+  else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
-+
+  else if(variable == "d0Sign"){
-+
+    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
-+
+    if(d0Sign < 0) value = -0.5;
-+
+    else if (d0Sign > 0) value = 0.5;
-+
+    else value = -999;
-+
+  }
-+
+  else if(variable == "chargeProduct"){
-+
+    value = object1->charge*object2->charge;
-+
+  }
-+
+  else if(variable == "electronCorrectedD0Vertex"){
-+
+    value = object1->correctedD0Vertex;
-+
+  }
-+
+  else if(variable == "muonCorrectedD0Vertex"){
-+
+    value = object2->correctedD0Vertex;
-+
+  }
-+
+  else if(variable == "electronCorrectedD0"){
-+
+    value = object1->correctedD0;
-+
+  }
-+
+  else if(variable == "muonCorrectedD0"){
-+
+    value = object2->correctedD0;
-+
+  }
-+
+  else if(variable == "electronDetIso"){
-+
+    value = (object1->trackIso) / object1->pt;
-+
+  }
-+
+  else if(variable == "muonDetIso"){
-+
+    value = (object2->trackIsoDR03) / object2->pt;
-+
+  }
-+
+  else if(variable == "chargeProduct"){
-+
+    value = object1->charge * object2->charge;
-+
+  }
-+
-+
-+
+  else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
-+
-+
+  value = applyFunction(function, value);
-+
-+
+  return value;
-+
+}
-+
-+
-+
+// Calculate the number of tracks in cone of DeltaR<0.5 around track1.
-+
+// Return true iff no other tracks are found in this cone.
-+
+int
-+
+OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){
-+
+  for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){
-+
+    if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself.
-+
+    double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi);
-+
+    if(deltaRtrk < 0.5) return 0;
-+
+  }
-+
+  return 1;
-+
-+
+}
-+
-+
-+
+double
-+
+OSUAnalysis::getTrkPtRes (const BNtrack* track1){
-+
-+
+  double ptTrue = getTrkPtTrue(track1, mcparticles.product());
-+
+  double PtRes = (track1->pt - ptTrue) / ptTrue;
-+
-+
+  return PtRes;
-+
-+
+}
-+
-+
-+
+double
-+
+OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){
-+
+  double value = -99;
-+
+  double genDeltaRLowest = 999;
-+
-+
+  for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){
-+
+    double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi);
-+
+    if (genDeltaRtemp < genDeltaRLowest) {
-+
+      genDeltaRLowest = genDeltaRtemp;
-+
+      if (genDeltaRLowest < 0.05) {   // Only consider it truth-matched if DeltaR<0.15.
-+
+        double ptTrue = genPart->pt;
-+
+        value = ptTrue;
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  return value;
-+
-+
+}
-+
-+
+//creates a map of the dead Ecal channels in the barrel and endcap
-+
+//to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here:
-+
+//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup
-+
+void
-+
+OSUAnalysis::WriteDeadEcal (){
-+
+  double etaEcal, phiEcal;
-+
+  ifstream DeadEcalFile(deadEcalFile_);
-+
+  if(!DeadEcalFile) {
-+
+    cout << "Error: DeadEcalFile has not been found." << endl;
-+
+    return;
-+
+  }
-+
+  if(DeadEcalVec.size()!= 0){
-+
+    cout << "Error: DeadEcalVec has a nonzero size" << endl;
-+
+    return;
-+
+  }
-+
+  while(!DeadEcalFile.eof())
-+
+    {
-+
+      DeadEcalFile >> etaEcal >> phiEcal;
-+
+      DeadEcal newChan;
-+
+      newChan.etaEcal = etaEcal;
-+
+      newChan.phiEcal = phiEcal;
-+
+      DeadEcalVec.push_back(newChan);
-+
+    }
-+
+  if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl;
-+
+}
-+
-+
+//if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0
-+
+int
-+
+OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){
-+
+  double deltaRLowest = 999;
-+
+  int value = 0;
-+
+  if (DeadEcalVec.size() == 0) WriteDeadEcal();
-+
+  for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){
-+
+    double eta = ecal->etaEcal;
-+
+    double phi = ecal->phiEcal;
-+
+    double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi);
-+
+    if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp;
-+
+  }
-+
+  if (deltaRLowest<0.05) {value = 1;}
-+
+  else {value = 0;}
-+
+  return value;
-+
+}
-+
-+
-+
-+
-+
+double
+OSUAnalysis::applyFunction(string function, double value){
+  if(function == "abs") value = fabs(value);
-+
+  else if(function == "fabs") value = fabs(value);
-+
+  else if(function == "log10") value = log10(value);
-+
+  else if(function == "log") value = log10(value);
-+
+  else if(function == "") value = value;
-+
+  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
+  return value;
+  for (uint object = 0; object != inputCollection->size(); object++){
-+
+    bool decision = true;//object passes if this cut doesn't cut on that type of object
-–
+    if(currentCut.inputCollection == inputType){
-<
+      double value = valueLookup(&inputCollection->at(object), currentCut.variable, currentCut.function);
->
+    if(currentCut.inputCollection == inputType){
-<
+      decision = evaluateComparison(value,currentCut.comparativeOperator,currentCut.cutValue);
->
+      vector<bool> subcutDecisions;
->
+      for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
->
+        double value = valueLookup(&inputCollection->at(object), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
->
+        subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
->
+      }
->
+      if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
->
+      else{
->
+        bool tempDecision = true;
->
+        for( int subcutIndex = 0;subcutIndex != currentCut.numSubcuts-1; subcutIndex++){
->
+          if(currentCut.logicalOperators.at(subcutIndex) == "&" || currentCut.logicalOperators.at(subcutIndex) == "&&")
->
+            tempDecision = subcutDecisions.at(subcutIndex) && subcutDecisions.at(subcutIndex+1);
->
+          else if(currentCut.logicalOperators.at(subcutIndex) == "|"|| currentCut.logicalOperators.at(subcutIndex) == "||")
->
+            tempDecision = subcutDecisions.at(subcutIndex) || subcutDecisions.at(subcutIndex+1);
->
+        }
->
+        decision = tempDecision;
->
+      }
-<
+    individualFlags[inputType].at(currentCutIndex).push_back(decision);
->
+    individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
+    //set flags for objects that pass each cut AND all the previous cuts
+    bool previousCumulativeFlag = true;
+    for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
-<
+      if(previousCumulativeFlag && individualFlags[inputType].at(previousCutIndex).at(object)) previousCumulativeFlag = true;
->
+      if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(object)) previousCumulativeFlag = true;
+      else{ previousCumulativeFlag = false; break;}
-<
+    cumulativeFlags[inputType].at(currentCutIndex).push_back(previousCumulativeFlag && decision);
->
+    cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
->
->
->
+  }
->
->
+}
->
->
->
+template <class InputCollection1, class InputCollection2>
->
+void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
->
+                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
->
->
->
+  bool sameObjects = false;
->
+  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
->
->
->
+  int counter = 0;
->
+  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
->
+    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
->
->
+      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
->
->
->
+      bool decision = true;//object passes if this cut doesn't cut on that type of object
->
->
+      if(currentCut.inputCollection == inputType){
->
->
+        vector<bool> subcutDecisions;
->
+        for( int subcutIndex = 0; subcutIndex != currentCut.numSubcuts; subcutIndex++){
->
+          double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), currentCut.variables.at(subcutIndex), currentCut.functions.at(subcutIndex));
->
+          subcutDecisions.push_back(evaluateComparison(value,currentCut.comparativeOperators.at(subcutIndex),currentCut.cutValues.at(subcutIndex)));
->
+        }
->
->
+        if(currentCut.numSubcuts == 1) decision = subcutDecisions.at(0);
->
+        else{
->
+          bool tempDecision = subcutDecisions.at(0);
->
+          for( int subcutIndex = 1; subcutIndex < currentCut.numSubcuts; subcutIndex++){
->
+            if(currentCut.logicalOperators.at(subcutIndex-1) == "&" || currentCut.logicalOperators.at(subcutIndex-1) == "&&")
->
+              tempDecision = tempDecision && subcutDecisions.at(subcutIndex);
->
+            else if(currentCut.logicalOperators.at(subcutIndex-1) == "|"|| currentCut.logicalOperators.at(subcutIndex-1) == "||")
->
+              tempDecision = tempDecision || subcutDecisions.at(subcutIndex);
->
+          }
->
+          decision = tempDecision;
->
+        }
->
+      }
->
+      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
->
->
+      //set flags for objects that pass each cut AND all the previous cuts
->
+      bool previousCumulativeFlag = true;
->
+      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
->
+        if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
->
+        else{ previousCumulativeFlag = false; break;}
->
+      }
->
+      //apply flags for the components of the composite object as well
->
+      bool currentCumulativeFlag = true;
->
+      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
->
+      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
->
+      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
->
+      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
->
+      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
->
->
+      counter++;
->
+    }
-+
+template <class InputCollection>
-<
+void OSUAnalysis::fillHistogram(TH1* histo, histogram parameters, InputCollection inputCollection){
->
+void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
->
->
+  for (uint object = 0; object != inputCollection->size(); object++){
->
->
+    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
->
->
+    string currentString = parameters.inputVariables.at(0);
->
+    string inputVariable = "";
->
+    string function = "";
->
+    if(currentString.find("(")==std::string::npos){
->
+      inputVariable = currentString;// variable to cut on
->
+    }
->
+    else{
->
+      function = currentString.substr(0,currentString.find("("));//function comes before the "("
->
+      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
->
+      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
->
+    }
-+
+    double value = valueLookup(&inputCollection->at(object), inputVariable, function);
-+
+    histo->Fill(value,scaleFactor);
-+
-+
+  }
-+
-+
+}
-+
-+
+template <class InputCollection1, class InputCollection2>
-+
+void OSUAnalysis::fill1DHistogram(TH1* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
-+
-+
+  bool sameObjects = false;
-+
+  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
-+
-+
+  int pairCounter = 0;
-+
+  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
-+
+    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
-+
-+
+      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
-+
-+
+      //only take objects which have passed all cuts and pairs which have passed all cuts
-+
+      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
-+
+      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
-+
+      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
-+
-+
+      string currentString = parameters.inputVariables.at(0);
-+
+      string inputVariable = "";
-+
+      string function = "";
-+
+      if(currentString.find("(")==std::string::npos){
-+
+        inputVariable = currentString;// variable to cut on
-+
+      }
-+
+      else{
-+
+        function = currentString.substr(0,currentString.find("("));//function comes before the "("
-+
+        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
-+
+        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
-+
+      }
-+
-+
+      double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
-+
+      histo->Fill(value,scaleFactor);
-+
-+
+      pairCounter++;
-+
+    }
-+
+  }
-+
-+
+}
-+
-+
-+
+template <class InputCollection>
-+
+void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection inputCollection,vector<bool> flags, double scaleFactor){
+  for (uint object = 0; object != inputCollection->size(); object++){
-<
+    if(parameters.inputVariables.size() == 1){
-<
+      double value = valueLookup(&inputCollection->at(object), parameters.inputVariables.at(0), parameters.function);
-<
+      histo->Fill(value);
-<
+    }
-<
+    else if(parameters.inputVariables.size() == 2){
-<
+      double valueX = valueLookup(&inputCollection->at(object), parameters.inputVariables.at(0), parameters.function);
-<
+      double valueY = valueLookup(&inputCollection->at(object), parameters.inputVariables.at(1), parameters.function);
-<
+      histo->Fill(valueX,valueY);
->
->
+    if(!plotAllObjectsInPassingEvents_ && !flags.at(object)) continue;
->
->
+    string currentString = parameters.inputVariables.at(0);
->
+    string inputVariable = "";
->
+    string function = "";
->
+    if(currentString.find("(")==std::string::npos){
->
+      inputVariable = currentString;// variable to cut on
-+
+    else{
-+
+      function = currentString.substr(0,currentString.find("("));//function comes before the "("
-+
+      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
-+
+      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
-+
+    }
-+
+    double valueX = valueLookup(&inputCollection->at(object), inputVariable, function);
-+
-+
+    currentString = parameters.inputVariables.at(1);
-+
+    inputVariable = "";
-+
+    function = "";
-+
+    if(currentString.find("(")==std::string::npos){
-+
+      inputVariable = currentString;// variable to cut on
-+
+    }
-+
+    else{
-+
+      function = currentString.substr(0,currentString.find("("));//function comes before the "("
-+
+      inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
-+
+      inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
-+
+    }
-+
-+
+    double valueY = valueLookup(&inputCollection->at(object), inputVariable, function);
-+
-+
+    histo->Fill(valueX,valueY,scaleFactor);
-+
-+
+template <class InputCollection1, class InputCollection2>
-+
+void OSUAnalysis::fill2DHistogram(TH2* histo, histogram parameters, InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, vector<bool> pairFlags, double scaleFactor){
-+
-+
+  bool sameObjects = false;
-+
+  if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
-+
-+
+  int pairCounter = 0;
-+
+  for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
-+
+    for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
-+
-+
+      if(sameObjects && object1 >= object2) continue;//account for duplicate pairs if both collections are the same
-+
-+
+      //only take objects which have passed all cuts and pairs which have passed all cuts
-+
+      if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
-+
+      if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
-+
+      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
-+
-+
+      string currentString = parameters.inputVariables.at(0);
-+
+      string inputVariable = "";
-+
+      string function = "";
-+
+      if(currentString.find("(")==std::string::npos){
-+
+        inputVariable = currentString;// variable to cut on
-+
+      }
-+
+      else{
-+
+        function = currentString.substr(0,currentString.find("("));//function comes before the "("
-+
+        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
-+
+        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
-+
+      }
-+
+      double valueX = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
-+
-+
+      currentString = parameters.inputVariables.at(1);
-+
+      inputVariable = "";
-+
+      function = "";
-+
+      if(currentString.find("(")==std::string::npos){
-+
+        inputVariable = currentString;// variable to cut on
-+
+      }
-+
+      else{
-+
+        function = currentString.substr(0,currentString.find("("));//function comes before the "("
-+
+        inputVariable = currentString.substr(currentString.find("(")+1);//get rest of string
-+
+        inputVariable = inputVariable.substr(0,inputVariable.size()-1);//remove trailing ")"
-+
+      }
-+
+      double valueY = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
-+
+      histo->Fill(valueX,valueY,scaleFactor);
-<
+DEFINE_FWK_MODULE(OSUAnalysis);
->
+      pairCounter++;
->
->
+    }
->
+  }
->
->
+}
->
->
->
+template <class InputObject>
->
+int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
-+
+  int bestMatchIndex = -1;
-+
+  double bestMatchDeltaR = 999;
-+
+  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
-+
-+
+    double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
-+
+    if(currentDeltaR > 0.05) continue;
-+
+//     cout << std::setprecision(3) << std::setw(20)
-+
+//          << "\tcurrentParticle:  eta = " << mcparticles->at(mcparticle - mcparticles->begin()).eta
-+
+//          << std::setw(20)
-+
+//          << "\tphi = " << mcparticles->at(mcparticle - mcparticles->begin()).phi
-+
+//          << std::setw(20)
-+
+//          << "\tdeltaR = " << currentDeltaR
-+
+//          << std::setprecision(1)
-+
+//          << std::setw(20)
-+
+//          << "\tid = " << mcparticles->at(mcparticle - mcparticles->begin()).id
-+
+//          << std::setw(20)
-+
+//          << "\tmotherId = " << mcparticles->at(mcparticle - mcparticles->begin()).motherId
-+
+//          << std::setw(20)
-+
+//          << "\tstatus = " << mcparticles->at(mcparticle - mcparticles->begin()).status<< endl;
-+
+    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
-+
+      bestMatchIndex = mcparticle - mcparticles->begin();
-+
+      bestMatchDeltaR = currentDeltaR;
-+
+    }
-+
-+
+  }
-+
+//   if(bestMatchDeltaR != 999)  cout << "bestMatch:  deltaR = " << bestMatchDeltaR << "   id = " << mcparticles->at(bestMatchIndex).id << "   motherId = " << mcparticles->at(bestMatchIndex).motherId << endl;
-+
+//   else cout << "no match found..." << endl;
-+
+  return bestMatchIndex;
-+
-+
+}
-+
-+
-+
+int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
-+
-+
+  int bestMatchIndex = -1;
-+
+  double bestMatchDeltaR = 999;
-+
-+
+  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
-+
-+
+    if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
-+
-+
+    double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
-+
+    if(currentDeltaR > 0.05) continue;
-+
-+
+    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
-+
+      bestMatchIndex = mcparticle - mcparticles->begin();
-+
+      bestMatchDeltaR = currentDeltaR;
-+
+    }
-+
-+
+  }
-+
-+
+  return bestMatchIndex;
-+
+}
-+
-+
-+
+// bin      particle type
-+
+// ---      -------------
-+
+//  0        unmatched
-+
+//  1        u
-+
+//  2        d
-+
+//  3        s
-+
+//  4        c
-+
+//  5        b
-+
+//  6        t
-+
+//  7        e
-+
+//  8        mu
-+
+//  9        tau
-+
+// 10        nu
-+
+// 11        g
-+
+// 12        gamma
-+
+// 13        Z
-+
+// 14        W
-+
+// 15        light meson
-+
+// 16        K meson
-+
+// 17        D meson
-+
+// 18        B meson
-+
+// 19        light baryon
-+
+// 20        strange baryon
-+
+// 21        charm baryon
-+
+// 22        bottom baryon
-+
+// 23        other
-+
-+
+int OSUAnalysis::getPdgIdBinValue(int pdgId){
-+
-+
+  int binValue = -999;
-+
+  int absPdgId = fabs(pdgId);
-+
+  if(pdgId == -1) binValue = 0;
-+
+  else if(absPdgId <= 6 ) binValue = absPdgId;
-+
+  else if(absPdgId == 11 ) binValue = 7;
-+
+  else if(absPdgId == 13 ) binValue = 8;
-+
+  else if(absPdgId == 15 ) binValue = 9;
-+
+  else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
-+
+  else if(absPdgId == 21 ) binValue = 11;
-+
+  else if(absPdgId == 22 ) binValue = 12;
-+
+  else if(absPdgId == 23 ) binValue = 13;
-+
+  else if(absPdgId == 24 ) binValue = 14;
-+
+  else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130  ) binValue = 15;
-+
+  else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400)  ) binValue = 16;
-+
+  else if(absPdgId > 400 && absPdgId < 500  ) binValue = 17;
-+
+  else if(absPdgId > 500 && absPdgId < 600  ) binValue = 18;
-+
+  else if(absPdgId > 1000 && absPdgId < 3000  ) binValue = 19;
-+
+  else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
-+
+  else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
-+
+  else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
-+
-+
+  else binValue = 23;
-+
-+
+  return binValue;
-+
-+
+}
-+
-+
+const BNprimaryvertex *
-+
+OSUAnalysis::chosenVertex ()
-+
+{
-+
+  const BNprimaryvertex *chosenVertex = 0;
-+
+  if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
-+
+    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back().size() ? cumulativeFlags.at("primaryvertexs").back() :
-+
+                               cumulativeFlags.at("primaryvertexs").at(cumulativeFlags.at("primaryvertexs").size() - 2);
-+
+    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
-+
+      if(!vertexFlags.at(vertexIndex)) continue;
-+
+      chosenVertex = & primaryvertexs->at(vertexIndex);
-+
+      break;
-+
+    }
-+
+  }
-+
+  else {
-+
+    chosenVertex = &primaryvertexs->at(0);
-+
+  }
-+
-+
+  return chosenVertex;
-+
+}
-+
-+
+const BNmet *
-+
+OSUAnalysis::chosenMET ()
-+
+{
-+
+  const BNmet *chosenMET = 0;
-+
+  if(std::find(objectsToCut.begin(), objectsToCut.end(), "mets") != objectsToCut.end()) {
-+
+    vector<bool> metFlags = cumulativeFlags.at("mets").back().size() ? cumulativeFlags.at("mets").back() :
-+
+                            cumulativeFlags.at("mets").at(cumulativeFlags.at("mets").size() - 2);
-+
+    for (uint metIndex = 0; metIndex != metFlags.size(); metIndex++){
-+
+      if(!metFlags.at(metIndex)) continue;
-+
+      chosenMET = & mets->at(metIndex);
-+
+      break;
-+
+    }
-+
+  }
-+
+  else {
-+
+    chosenMET = &mets->at(0);
-+
+  }
-+
-+
+  return chosenMET;
-+
+}
-+
-+
+const BNelectron *
-+
+OSUAnalysis::chosenElectron ()
-+
+{
-+
+  const BNelectron *chosenElectron = 0;
-+
+  if(std::find(objectsToCut.begin(), objectsToCut.end(), "electrons") != objectsToCut.end()) {
-+
+    vector<bool> electronFlags = cumulativeFlags.at("electrons").back().size() ? cumulativeFlags.at("electrons").back() :
-+
+                                 cumulativeFlags.at("electrons").at(cumulativeFlags.at("electrons").size() - 2);
-+
+    for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
-+
+      if(!electronFlags.at(electronIndex)) continue;
-+
+      chosenElectron = & electrons->at(electronIndex);
-+
+      break;
-+
+    }
-+
+  }
-+
+  else {
-+
+    chosenElectron = &electrons->at(0);
-+
+  }
-+
-+
+  return chosenElectron;
-+
+}
-+
-+
+const BNmuon *
-+
+OSUAnalysis::chosenMuon ()
-+
+{
-+
+  const BNmuon *chosenMuon = 0;
-+
+  if(std::find(objectsToCut.begin(), objectsToCut.end(), "muons") != objectsToCut.end()) {
-+
+    vector<bool> muonFlags = cumulativeFlags.at("muons").back().size() ? cumulativeFlags.at("muons").back() :
-+
+                             cumulativeFlags.at("muons").at(cumulativeFlags.at("muons").size() - 2);
-+
+    for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
-+
+      if(!muonFlags.at(muonIndex)) continue;
-+
+      chosenMuon = & muons->at(muonIndex);
-+
+      break;
-+
+    }
-+
+  }
-+
+  else {
-+
+    chosenMuon = &muons->at(0);
-+
+  }
-+
-+
+  return chosenMuon;
-+
+}
-+
-+
-+
+DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc (file contents): Revision 1.12 by ahart, Fri Feb 8 11:07:51 2013 UTC vs. Revision 1.32 by qpython, Tue Mar 19 18:05:23 2013 UTC

Diff Legend

Comparing UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc (file contents):
Revision 1.12 by ahart, Fri Feb 8 11:07:51 2013 UTC vs.
Revision 1.32 by qpython, Tue Mar 19 18:05:23 2013 UTC