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

Comparing UserCode/OSUT3Analysis/AnaTools/plugins/OSUAnalysis.cc (file contents):
Revision 1.18 by lantonel, Fri Feb 22 10:38:08 2013 UTC vs.
Revision 1.23 by lantonel, Mon Mar 4 22:38:39 2013 UTC

# Line 37 | Line 37 | OSUAnalysis::OSUAnalysis (const edm::Par
37    std::vector<TFileDirectory> directories;
38  
39    //always get vertex collection so we can assign the primary vertex in the event
40 <  allNecessaryObjects.push_back("primaryvertexs");
40 >  objectsToGet.push_back("primaryvertexs");
41    //always make the plot of number of primary verticex (to check pile-up reweighting)
42    objectsToPlot.push_back("primaryvertexs");
43 +  //always make flags for the vertices to pick the best good one
44 +  //  objectsToCut.push_back("primaryvertexs");
45 +
46 +  //always get the MC particles to do GEN-matching
47 +  objectsToGet.push_back("mcparticles");
48  
49    //always get the event collection to do pile-up reweighting
50 <  allNecessaryObjects.push_back("events");
50 >  objectsToGet.push_back("events");
51  
52    //parse the histogram definitions
53    for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
# Line 50 | Line 55 | OSUAnalysis::OSUAnalysis (const edm::Par
55      string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
56      if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
57      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
58 +      objectsToGet.push_back(tempInputCollection);
59        objectsToPlot.push_back(tempInputCollection);
54      allNecessaryObjects.push_back(tempInputCollection);
60      }
61 <    else{//pair of objects, need to add them both to the things to allNecessaryObjects
61 >    else{//pair of objects, need to add them both to the things to objectsToGet
62        int dashIndex = tempInputCollection.find("-");
63        int spaceIndex = tempInputCollection.find(" ");
64        int secondWordLength = spaceIndex - dashIndex;
65 <      allNecessaryObjects.push_back(tempInputCollection);
66 <      allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s");
67 <      allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
65 >      objectsToGet.push_back(tempInputCollection);
66 >      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
67 >      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
68        objectsToPlot.push_back(tempInputCollection);
69        objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
70        objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
71 +      objectsToCut.push_back(tempInputCollection);
72 +      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
73 +      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
74 +
75        }
76  
77      vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
# Line 80 | Line 89 | OSUAnalysis::OSUAnalysis (const edm::Par
89  
90      }
91    }
92 <  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object
92 >  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
93    sort( objectsToPlot.begin(), objectsToPlot.end() );
94    objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
95  
96  
97  
98 +  //add histograms with the gen-matched id, mother id, and grandmother id
99 +  for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
100 +
101 +    string currentObject = objectsToPlot.at(currentObjectIndex);
102 +    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
103 +
104 +    histogram tempIdHisto;
105 +    histogram tempMomIdHisto;
106 +    histogram tempGmaIdHisto;
107 +
108 +    tempIdHisto.inputCollection = currentObject;
109 +    tempMomIdHisto.inputCollection = currentObject;
110 +    tempGmaIdHisto.inputCollection = currentObject;
111 +
112 +    currentObject = currentObject.substr(0, currentObject.size()-1);
113 +    tempIdHisto.name = currentObject+"GenMatchId";
114 +    tempMomIdHisto.name = currentObject+"GenMatchMotherId";
115 +    tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
116 +
117 +    currentObject.at(0) = toupper(currentObject.at(0));
118 +    tempIdHisto.title = currentObject+" Gen-matched Particle";
119 +    tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
120 +    tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
121 +
122 +    int maxNum = 24;
123 +    vector<double> binVector;
124 +    binVector.push_back(maxNum);
125 +    binVector.push_back(0);
126 +    binVector.push_back(maxNum);
127 +
128 +    tempIdHisto.bins = binVector;
129 +    tempIdHisto.inputVariables.push_back("genMatchedId");
130 +    tempMomIdHisto.bins = binVector;
131 +    tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
132 +    tempGmaIdHisto.bins = binVector;
133 +    tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
134 +
135 +    histograms.push_back(tempIdHisto);
136 +    histograms.push_back(tempMomIdHisto);
137 +    histograms.push_back(tempGmaIdHisto);
138 +    
139 +  }
140 +
141 +
142  
143    channel tempChannel;
144    //loop over all channels (event selections)
# Line 104 | Line 157 | OSUAnalysis::OSUAnalysis (const edm::Par
157        triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
158        for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
159          tempChannel.triggers.push_back(triggerNames.at(trigger));
160 <      allNecessaryObjects.push_back("triggers");
160 >      objectsToGet.push_back("triggers");
161      }
162  
163  
# Line 146 | Line 199 | OSUAnalysis::OSUAnalysis (const edm::Par
199        }
200  
201  
202 +      if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
203 +
204 + // bin      particle type
205 + // ---      -------------
206 + //  0        unmatched
207 + //  1        u
208 + //  2        d
209 + //  3        s
210 + //  4        c
211 + //  5        b
212 + //  6        t
213 + //  7        e
214 + //  8        mu
215 + //  9        tau
216 + // 10        nu
217 + // 11        g
218 + // 12        gamma
219 + // 13        Z
220 + // 14        W
221 + // 15        light meson
222 + // 16        K meson
223 + // 17        D meson
224 + // 18        B meson
225 + // 19        light baryon
226 + // 20        strange baryon
227 + // 21        charm baryon
228 + // 22        bottom baryon
229 + // 23        other
230 +
231 +      vector<TString> labelArray;
232 +      labelArray.push_back("unmatched");
233 +      labelArray.push_back("u");
234 +      labelArray.push_back("d");
235 +      labelArray.push_back("s");
236 +      labelArray.push_back("c");
237 +      labelArray.push_back("b");
238 +      labelArray.push_back("t");
239 +      labelArray.push_back("e");
240 +      labelArray.push_back("#mu");
241 +      labelArray.push_back("#tau");
242 +      labelArray.push_back("#nu");
243 +      labelArray.push_back("g");
244 +      labelArray.push_back("#gamma");
245 +      labelArray.push_back("Z");
246 +      labelArray.push_back("W");
247 +      labelArray.push_back("light meson");
248 +      labelArray.push_back("K meson");
249 +      labelArray.push_back("D meson");
250 +      labelArray.push_back("B meson");
251 +      labelArray.push_back("light baryon");
252 +      labelArray.push_back("strange baryon");
253 +      labelArray.push_back("charm baryon");
254 +      labelArray.push_back("bottom baryon");
255 +      labelArray.push_back("other");
256 +
257 +      for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
258 +        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
259 +      }
260      }
261      //book a histogram for the number of each object type to be plotted
262      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
# Line 171 | Line 282 | OSUAnalysis::OSUAnalysis (const edm::Par
282      }
283  
284  
285 +    
286 +
287 +
288      //get list of cuts for this channel
289      vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
290  
# Line 181 | Line 295 | OSUAnalysis::OSUAnalysis (const edm::Par
295        string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
296        tempCut.inputCollection = tempInputCollection;
297        if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
298 <        allNecessaryObjects.push_back(tempInputCollection);
298 >        objectsToGet.push_back(tempInputCollection);
299 >        objectsToCut.push_back(tempInputCollection);
300        }
301 <      else{//pair of objects, need to add them both to the things to allNecessaryObjects
301 >      else{//pair of objects, need to add them both to the things to objectsToGet
302          int dashIndex = tempInputCollection.find("-");
303          int spaceIndex = tempInputCollection.find(" ");
304          int secondWordLength = spaceIndex - dashIndex;
305 <        allNecessaryObjects.push_back(tempInputCollection);
306 <        allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s");
307 <        allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
305 >        objectsToGet.push_back(tempInputCollection);
306 >        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
307 >        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
308 >        objectsToCut.push_back(tempInputCollection);
309 >        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
310 >        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
311  
312        }
313  
# Line 259 | Line 377 | OSUAnalysis::OSUAnalysis (const edm::Par
377    }//end loop over channels
378  
379  
380 <  //make unique vector of objects we need to cut on (so we make sure to get them from the event)
381 <  sort( allNecessaryObjects.begin(), allNecessaryObjects.end() );
382 <  allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() );
380 >  //make unique vector of objects we need to get from the event
381 >  sort( objectsToGet.begin(), objectsToGet.end() );
382 >  objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
383 >  //make unique vector of objects we need to cut on
384 >  sort( objectsToCut.begin(), objectsToCut.end() );
385 >  objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
386  
387  
388   }
# Line 279 | Line 400 | void
400   OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
401   {
402  
282
403    // Retrieve necessary collections from the event.
404 <  edm::Handle<BNtriggerCollection> triggers;
405 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "triggers") != allNecessaryObjects.end())
404 >
405 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
406      event.getByLabel (triggers_, triggers);
407  
408 <  edm::Handle<BNjetCollection> jets;
289 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end())
408 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
409      event.getByLabel (jets_, jets);
410  
411 <  edm::Handle<BNmuonCollection> muons;
293 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end())
411 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
412      event.getByLabel (muons_, muons);
413  
414 <  edm::Handle<BNelectronCollection> electrons;
297 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end())
414 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
415      event.getByLabel (electrons_, electrons);
416  
417 <  edm::Handle<BNeventCollection> events;
301 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end())
417 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
418      event.getByLabel (events_, events);
419  
420 <  edm::Handle<BNtauCollection> taus;
305 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end())
420 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
421      event.getByLabel (taus_, taus);
422  
423 <  edm::Handle<BNmetCollection> mets;
309 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end())
423 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
424      event.getByLabel (mets_, mets);
425  
426 <  edm::Handle<BNtrackCollection> tracks;
313 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end())
426 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
427      event.getByLabel (tracks_, tracks);
428  
429 <  edm::Handle<BNgenjetCollection> genjets;
317 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end())
429 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
430      event.getByLabel (genjets_, genjets);
431  
432 <  edm::Handle<BNmcparticleCollection> mcparticles;
321 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end())
432 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
433      event.getByLabel (mcparticles_, mcparticles);
434  
435 <  edm::Handle<BNprimaryvertexCollection> primaryvertexs;
325 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end())
435 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
436      event.getByLabel (primaryvertexs_, primaryvertexs);
437  
438 <  edm::Handle<BNbxlumiCollection> bxlumis;
329 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end())
438 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
439      event.getByLabel (bxlumis_, bxlumis);
440  
441 <  edm::Handle<BNphotonCollection> photons;
333 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end())
441 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
442      event.getByLabel (photons_, photons);
443  
444 <  edm::Handle<BNsuperclusterCollection> superclusters;
337 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end())
444 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
445      event.getByLabel (superclusters_, superclusters);
446  
447  
448    //get pile-up event weight
449 <  double puScaleFactor = 0;
449 >  double puScaleFactor = 1.00;
450    if(datasetType_ != "data")
451      puScaleFactor = puWeight_->at (events->at (0).numTruePV);
345  else
346    puScaleFactor = 1.00;
347
452  
453    //loop over all channels
454  
# Line 362 | Line 466 | OSUAnalysis::analyze (const edm::Event &
466      }
467  
468      //loop over all cuts
469 +
470 +
471      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
472        cut currentCut = currentChannel.cuts.at(currentCutIndex);
473  
474 <      for(uint currentObjectIndex = 0; currentObjectIndex != allNecessaryObjects.size(); currentObjectIndex++){
475 <        string currentObject = allNecessaryObjects.at(currentObjectIndex);
474 >      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
475 >        string currentObject = objectsToCut.at(currentObjectIndex);
476  
477 +        //initialize maps to get ready to set cuts
478          individualFlags[currentObject].push_back (vector<bool> ());
479          cumulativeFlags[currentObject].push_back (vector<bool> ());
480  
481 +      }
482 +      for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){
483 +        string currentObject = objectsToCut.at(currentObjectIndex);
484 +
485 +        int flagsForPairCutsIndex = max(int(currentCutIndex-1),0);
486  
487          if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets");
488          else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons");
# Line 386 | Line 498 | OSUAnalysis::analyze (const edm::Event &
498          else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons");
499          else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters");
500  
501 <        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(),"muon-muon pairs");
502 <        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(),"electron-electron pairs");
503 <        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(),"electron-muon pairs");
501 >        
502 >        else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \
503 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
504 >                                                                   cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
505 >                                                                   "muon-muon pairs");
506 >        else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \
507 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
508 >                                                                           cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
509 >                                                                           "electron-electron pairs");
510 >        else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \
511 >                                                                       cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \
512 >                                                                       cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \
513 >                                                                       "electron-muon pairs");
514  
515        }
516  
# Line 397 | Line 519 | OSUAnalysis::analyze (const edm::Event &
519      }//end loop over all cuts
520  
521  
400
522      //use cumulative flags to apply cuts at event level
523  
524      bool eventPassedAllCuts = true;
# Line 430 | Line 551 | OSUAnalysis::analyze (const edm::Event &
551      if(!eventPassedAllCuts)continue;
552  
553  
554 <
554 >    
555  
556      //set position of primary vertex in event, in order to calculate quantities relative to it
557 <    primaryVertex_ = 0;
558 <    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
559 <    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
560 <      if(!vertexFlags.at(vertexIndex)) continue;
561 <      primaryVertex_ = new BNprimaryvertex (primaryvertexs->at (vertexIndex));
562 <      break;
557 >    if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
558 >      vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
559 >      for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
560 >        if(!vertexFlags.at(vertexIndex)) continue;
561 >        chosenPrimaryVertex = & primaryvertexs->at(vertexIndex);
562 >        break;
563 >      }
564 >    }
565 >    else {
566 >      chosenPrimaryVertex = & primaryvertexs->at(0);
567      }
568  
569  
445
570      //filling histograms
571      for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){
572        histogram currentHistogram = histograms.at(histogramIndex);
# Line 450 | Line 574 | OSUAnalysis::analyze (const edm::Event &
574        if(currentHistogram.inputVariables.size() == 1){
575          TH1D* histo;
576          histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
577 +
578          if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
579 <        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
579 >        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
580          else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),\
581                                                                                         cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(),\
582                                                                                         cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
# Line 501 | Line 626 | OSUAnalysis::analyze (const edm::Event &
626        }
627      }
628  
504
629      //fills histograms with the sizes of collections
630      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
631        string currentObject = objectsToPlot.at(currentObjectIndex);
632  
633 <      if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
634 <      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
635 <      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
636 <      string tempCurrentObject = currentObject;
633 >      string objectToPlot = "";
634 >
635 >      if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
636 >      else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
637 >      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
638 >      else objectToPlot = currentObject;
639 >      string tempCurrentObject = objectToPlot;
640        tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
641        string histoName = "num" + tempCurrentObject;
642  
643  
644 <      if(currentObject == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
645 <      else if(currentObject == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
646 <      else if(currentObject == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
647 <      else if(currentObject == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
648 <      else if(currentObject == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
649 <      else if(currentObject == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
650 <      else if(currentObject == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
651 <      else if(currentObject == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
652 <      else if(currentObject == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
653 <      else if(currentObject == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
654 <      else if(currentObject == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
655 <      else if(currentObject == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
656 <      else if(currentObject == "primaryvertexs"){
644 >
645 >
646 >      //set position of primary vertex in event, in order to calculate quantities relative to it
647 >      if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
648 >        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
649 >        int numToPlot = 0;
650 >        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
651 >          if(lastCutFlags.at(currentFlag)) numToPlot++;
652 >        }
653 >        if(objectToPlot == "primaryvertexs"){
654 >          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
655 >          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
656 >        }
657 >        else
658 >          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,puScaleFactor);        
659 >      }
660 >      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
661 >      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
662 >      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
663 >      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
664 >      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
665 >      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
666 >      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
667 >      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
668 >      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
669 >      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
670 >      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
671 >      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
672 >      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor);
673 >      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor);
674 >      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor);
675 >      else if(objectToPlot == "primaryvertexs"){
676          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
677          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
678        }
533      else if(currentObject == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor);
534      else if(currentObject == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor);
535      else if(currentObject == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor);
679  
680      }
681 +
682 +
683 +
684 +
685    } //end loop over channel
686  
687    masterCutFlow_->fillCutFlow(puScaleFactor);
# Line 736 | Line 883 | OSUAnalysis::valueLookup (const BNmuon*
883    else if(variable == "ecalIso") value = object->ecalIso;
884    else if(variable == "hcalIso") value = object->hcalIso;
885    else if(variable == "caloIso") value = object->caloIso;
886 <  else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
886 >  else if(variable == "trackIsDR03") value = object->trackIsoDR03;
887    else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
888    else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
889    else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
# Line 875 | Line 1022 | OSUAnalysis::valueLookup (const BNmuon*
1022    else if(variable == "time_ndof") value = object->time_ndof;
1023  
1024    //user-defined variables
1025 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1026 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1027 <  else if(variable == "detIso") value = (object->trackIso) / object->pt;
1025 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1026 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1027 >  else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1028    else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1029    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1030    else if(variable == "tightID") {
# Line 888 | Line 1035 | OSUAnalysis::valueLookup (const BNmuon*
1035        && object->numberOfMatchedStations > 1    \
1036        && fabs(object->correctedD0Vertex) < 0.2  \
1037        && fabs(object->correctedDZ) < 0.5        \
1038 <      && object->numberOfValidPixelHits > 0     \
1038 >      && object->numberOfValidPixelHits > 0             \
1039        && object->numberOfLayersWithMeasurement > 5;
893    
894
1040    }
1041    else if(variable == "tightIDdisplaced"){
1042      value = object->isGlobalMuon > 0            \
898      && object->isPFMuon > 0                   \
1043        && object->normalizedChi2 < 10            \
1044        && object->numberOfValidMuonHits > 0      \
1045        && object->numberOfMatchedStations > 1    \
# Line 903 | Line 1047 | OSUAnalysis::valueLookup (const BNmuon*
1047        && object->numberOfLayersWithMeasurement > 5;
1048    }
1049  
1050 +  else if(variable == "genMatchedId"){
1051 +    int index = getGenMatchedParticleIndex(object);
1052 +    if(index == -1) value = 0;
1053 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1054 +  }
1055 +  else if(variable == "genMatchedMotherId"){
1056 +    int index = getGenMatchedParticleIndex(object);
1057 +    if(index == -1) value = 0;
1058 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1059 +  }
1060 +  else if(variable == "genMatchedGrandmotherId"){
1061 +    int index = getGenMatchedParticleIndex(object);
1062 +    if(index == -1) value = 0;
1063 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1064 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1065 +      if(motherIndex == -1) value = 0;
1066 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1067 +    }
1068 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1069 +  }
1070 +
1071 +
1072 +
1073    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1074  
1075    value = applyFunction(function, value);
# Line 1073 | Line 1240 | OSUAnalysis::valueLookup (const BNelectr
1240    else if(variable == "passConvVeto") value = object->passConvVeto;
1241  
1242    //user-defined variables
1243 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1244 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1243 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1244 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1245    else if(variable == "detIso") value = (object->trackIso) / object->pt;
1246    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1247 +  else if(variable == "correctedD0VertexInEB"){
1248 +    if(fabs(object->eta) < 0.8) value = object->correctedD0Vertex;
1249 +    else value = -999;
1250 +  }
1251 +  else if(variable == "correctedD0VertexOutEB"){
1252 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0Vertex;
1253 +    else value = -999;
1254 +  }
1255 +  else if(variable == "correctedD0VertexEE"){
1256 +    if(object->isEE) value = object->correctedD0Vertex;
1257 +    else value = -999;
1258 +  }
1259 +
1260 +  else if(variable == "correctedD0BeamspotInEB"){
1261 +    if(fabs(object->eta) < 0.8) value = object->correctedD0;
1262 +    else value = -999;
1263 +  }
1264 +  else if(variable == "correctedD0BeamspotOutEB"){
1265 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0;
1266 +    else value = -999;
1267 +  }
1268 +  else if(variable == "correctedD0BeamspotEE"){
1269 +    if(object->isEE) value = object->correctedD0;
1270 +    else value = -999;
1271 +  }
1272 +
1273 +  else if(variable == "correctedD0OriginInEB"){
1274 +    if(fabs(object->eta) < 0.8) value = object->tkD0;
1275 +    else value = -999;
1276 +  }
1277 +  else if(variable == "correctedD0OriginOutEB"){
1278 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->tkD0;
1279 +    else value = -999;
1280 +  }
1281 +  else if(variable == "correctedD0OriginEE"){
1282 +    if(object->isEE) value = object->tkD0;
1283 +    else value = -999;
1284 +  }
1285 +
1286 +  else if(variable == "tightIDdisplaced"){
1287 +    if (object->isEB)
1288 +      {
1289 +        value = fabs(object->delEtaIn) < 0.004 \
1290 +          && fabs (object->delPhiIn) < 0.03 \
1291 +          && object->scSigmaIEtaIEta < 0.01 \
1292 +          && object->hadOverEm < 0.12 \
1293 +          && object->absInvEMinusInvPin < 0.05;
1294 +      }
1295 +    else
1296 +      {
1297 +        value = fabs (object->delEtaIn) < 0.005 \
1298 +          && fabs (object->delPhiIn) < 0.02 \
1299 +          && object->scSigmaIEtaIEta < 0.03 \
1300 +          && object->hadOverEm < 0.10 \
1301 +          && object->absInvEMinusInvPin < 0.05;
1302 +      }
1303 +  }
1304 +
1305 +  else if(variable == "genMatchedId"){
1306 +    int index = getGenMatchedParticleIndex(object);
1307 +    if(index == -1) value = 0;
1308 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1309 +  }
1310 +  else if(variable == "genMatchedMotherId"){
1311 +    int index = getGenMatchedParticleIndex(object);
1312 +    if(index == -1) value = 0;
1313 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1314 +  }
1315 +  else if(variable == "genMatchedGrandmotherId"){
1316 +    int index = getGenMatchedParticleIndex(object);
1317 +    if(index == -1) value = 0;
1318 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1319 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1320 +      if(motherIndex == -1) value = 0;
1321 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1322 +    }
1323 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1324 +  }
1325 +
1326  
1327  
1328    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
# Line 1210 | Line 1456 | OSUAnalysis::valueLookup (const BNtau* o
1456    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1457  
1458  
1459 +  else if(variable == "genMatchedId"){
1460 +    int index = getGenMatchedParticleIndex(object);
1461 +    if(index == -1) value = 0;
1462 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1463 +  }
1464 +  else if(variable == "genMatchedMotherId"){
1465 +    int index = getGenMatchedParticleIndex(object);
1466 +    if(index == -1) value = 0;
1467 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1468 +  }
1469 +  else if(variable == "genMatchedGrandmotherId"){
1470 +    int index = getGenMatchedParticleIndex(object);
1471 +    if(index == -1) value = 0;
1472 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1473 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1474 +      if(motherIndex == -1) value = 0;
1475 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1476 +    }
1477 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1478 +  }
1479 +
1480 +
1481    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1482  
1483    value = applyFunction(function, value);
# Line 1313 | Line 1581 | OSUAnalysis::valueLookup (const BNtrack*
1581    else if(variable == "isHighPurity") value = object->isHighPurity;
1582  
1583  
1584 +  else if(variable == "genMatchedId"){
1585 +    int index = getGenMatchedParticleIndex(object);
1586 +    if(index == -1) value = 0;
1587 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1588 +  }
1589 +  else if(variable == "genMatchedMotherId"){
1590 +    int index = getGenMatchedParticleIndex(object);
1591 +    if(index == -1) value = 0;
1592 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1593 +  }
1594 +  else if(variable == "genMatchedGrandmotherId"){
1595 +    int index = getGenMatchedParticleIndex(object);
1596 +    if(index == -1) value = 0;
1597 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1598 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1599 +      if(motherIndex == -1) value = 0;
1600 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1601 +    }
1602 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1603 +  }
1604 +
1605 +
1606 +
1607    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1608  
1609    value = applyFunction(function, value);
# Line 1440 | Line 1731 | OSUAnalysis::valueLookup (const BNmcpart
1731  
1732    //user-defined variables
1733    else if (variable == "d0"){
1734 <    double vx = object->vx - primaryVertex_->x,
1735 <      vy = object->vy - primaryVertex_->y,
1734 >    double vx = object->vx - chosenPrimaryVertex->x,
1735 >      vy = object->vy - chosenPrimaryVertex->y,
1736        px = object->px,
1737        py = object->py,
1738        pt = object->pt;
1739      value = (-vx * py + vy * px) / pt;
1740    }
1741    else if (variable == "dz"){
1742 <    double vx = object->vx - primaryVertex_->x,
1743 <      vy = object->vy - primaryVertex_->y,
1744 <      vz = object->vz - primaryVertex_->z,
1742 >    double vx = object->vx - chosenPrimaryVertex->x,
1743 >      vy = object->vy - chosenPrimaryVertex->y,
1744 >      vz = object->vz - chosenPrimaryVertex->z,
1745        px = object->px,
1746        py = object->py,
1747        pz = object->pz,
1748        pt = object->pt;
1749      value = vz - (vx * px + vy * py)/pt * (pz/pt);
1750    }
1751 <
1751 >  else if(variable == "v0"){
1752 >    value = sqrt(object->vx*object->vx + object->vy*object->vy);
1753 >  }
1754 >  else if(variable == "deltaV0"){
1755 >    value = sqrt((object->vx-chosenPrimaryVertex->x)*(object->vx-chosenPrimaryVertex->x) + (object->vy-chosenPrimaryVertex->y)*(object->vy-chosenPrimaryVertex->y));
1756 >  }
1757 >  else if (variable == "deltaVx"){
1758 >    value = object->vx - chosenPrimaryVertex->x;
1759 >  }
1760 >  else if (variable == "deltaVy"){
1761 >    value = object->vy - chosenPrimaryVertex->y;
1762 >  }
1763 >  else if (variable == "deltaVz"){
1764 >    value = object->vz - chosenPrimaryVertex->z;
1765 >  }
1766  
1767  
1768    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
# Line 1587 | Line 1892 | OSUAnalysis::valueLookup (const BNphoton
1892    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
1893  
1894  
1895 +  else if(variable == "genMatchedId"){
1896 +    int index = getGenMatchedParticleIndex(object);
1897 +    if(index == -1) value = 0;
1898 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1899 +  }
1900 +  else if(variable == "genMatchedMotherId"){
1901 +    int index = getGenMatchedParticleIndex(object);
1902 +    if(index == -1) value = 0;
1903 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1904 +  }
1905 +  else if(variable == "genMatchedGrandmotherId"){
1906 +    int index = getGenMatchedParticleIndex(object);
1907 +    if(index == -1) value = 0;
1908 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1909 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1910 +      if(motherIndex == -1) value = 0;
1911 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1912 +    }
1913 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1914 +  }
1915 +
1916 +
1917    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1918  
1919    value = applyFunction(function, value);
# Line 1622 | Line 1949 | OSUAnalysis::valueLookup (const BNmuon*
1949  
1950    double value = 0.0;
1951  
1952 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
1952 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
1953 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
1954    else if(variable == "invMass"){
1955      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
1956      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
1957 <    value = (fourVector1 + fourVector2).M();
1958 <  }
1957 >
1958 >    value = (fourVector1 + fourVector2).M();}
1959 >  else if(variable == "threeDAngle")
1960 >    {
1961 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
1962 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
1963 >      value = (threeVector1.Angle(threeVector2));
1964 >    }
1965 >
1966 >
1967 >
1968    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
1969    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
1970 <  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
1971 <
1970 >  else if(variable == "d0Sign"){
1971 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
1972 >    if(d0Sign < 0) value = -0.5;
1973 >    else if (d0Sign > 0) value = 0.5;
1974 >    else value = -999;
1975 >  }
1976 >  else if(variable == "muon1CorrectedD0Vertex"){
1977 >    value = object1->correctedD0Vertex;
1978 >  }
1979 >  else if(variable == "muon2CorrectedD0Vertex"){
1980 >    value = object2->correctedD0Vertex;
1981 >  }
1982    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1983  
1984 +
1985 +
1986 +
1987 +
1988 +
1989 +
1990 +
1991    value = applyFunction(function, value);
1992  
1993    return value;
# Line 1644 | Line 1998 | OSUAnalysis::valueLookup (const BNelectr
1998  
1999    double value = 0.0;
2000  
2001 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
2001 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2002 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2003    else if(variable == "invMass"){
2004      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2005      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2006 <    value = (fourVector1 + fourVector2).M();
2007 <  }
2006 >
2007 >    value = (fourVector1 + fourVector2).M();}
2008 > else if(variable == "threeDAngle")
2009 >    {
2010 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2011 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2012 >      value = (threeVector1.Angle(threeVector2));
2013 >    }
2014 >  
2015 >
2016 >
2017 >
2018    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2019    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2020 +
2021    else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2022  
2023 +
2024 +
2025 +  else if(variable == "d0Sign"){
2026 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2027 +    if(d0Sign < 0) value = -0.5;
2028 +    else if (d0Sign > 0) value = 0.5;
2029 +    else value = -999;
2030 +  }
2031 +  else if(variable == "electron1CorrectedD0Vertex"){
2032 +    value = object1->correctedD0Vertex;
2033 +  }
2034 +  else if(variable == "electron2CorrectedD0Vertex"){
2035 +    value = object2->correctedD0Vertex;
2036 +  }
2037 +
2038    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2039  
2040    value = applyFunction(function, value);
# Line 1666 | Line 2047 | OSUAnalysis::valueLookup (const BNelectr
2047  
2048    double value = 0.0;
2049  
2050 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
2050 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2051 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2052    else if(variable == "invMass"){
2053      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2054      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2055 <    value = (fourVector1 + fourVector2).M();
2056 <  }
2055 >
2056 >    value = (fourVector1 + fourVector2).M();}
2057 > else if(variable == "threeDAngle")
2058 >    {
2059 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2060 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2061 >      value = (threeVector1.Angle(threeVector2));
2062 >    }
2063 >
2064 >
2065    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2066    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2067 <  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2067 >  else if(variable == "d0Sign"){
2068 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2069 >    if(d0Sign < 0) value = -0.5;
2070 >    else if (d0Sign > 0) value = 0.5;
2071 >    else value = -999;
2072 >  }
2073 >  else if(variable == "electronCorrectedD0Vertex"){
2074 >    value = object1->correctedD0Vertex;
2075 >  }
2076 >  else if(variable == "muonCorrectedD0Vertex"){
2077 >    value = object2->correctedD0Vertex;
2078 >  }
2079 >  else if(variable == "electronCorrectedD0"){
2080 >    value = object1->correctedD0;
2081 >  }
2082 >  else if(variable == "muonCorrectedD0"){
2083 >    value = object2->correctedD0;
2084 >  }
2085 >
2086  
2087    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2088  
# Line 1688 | Line 2096 | double
2096   OSUAnalysis::applyFunction(string function, double value){
2097  
2098    if(function == "abs") value = fabs(value);
2099 +  else if(function == "fabs") value = fabs(value);
2100 +
2101  
2102 +  else if(function == "") value = value;
2103 +  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
2104  
2105    return value;
2106  
# Line 1698 | Line 2110 | OSUAnalysis::applyFunction(string functi
2110   template <class InputCollection>
2111   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2112  
2113 +
2114    for (uint object = 0; object != inputCollection->size(); object++){
2115  
2116 +
2117      bool decision = true;//object passes if this cut doesn't cut on that type of object
2118  
2119 +
2120      if(currentCut.inputCollection == inputType){
2121  
2122        vector<bool> subcutDecisions;
# Line 1723 | Line 2138 | void OSUAnalysis::setObjectFlags(cut &cu
2138      }
2139      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2140  
2141 +
2142      //set flags for objects that pass each cut AND all the previous cuts
2143      bool previousCumulativeFlag = true;
2144      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
# Line 1731 | Line 2147 | void OSUAnalysis::setObjectFlags(cut &cu
2147      }
2148      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2149  
2150 +
2151    }
2152  
2153   }
# Line 1738 | Line 2155 | void OSUAnalysis::setObjectFlags(cut &cu
2155  
2156   template <class InputCollection1, class InputCollection2>
2157   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
2158 <                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, string inputType){
2158 >                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
2159  
2160  
2161    bool sameObjects = false;
2162    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2163  
2164 +
2165    int counter = 0;  
2166    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2167      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
# Line 1781 | Line 2199 | void OSUAnalysis::setObjectFlags(cut &cu
2199          if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
2200          else{ previousCumulativeFlag = false; break;}
2201        }
2202 <      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2202 >      //apply flags for the components of the composite object as well
2203 >      bool currentCumulativeFlag = true;
2204 >      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
2205 >      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
2206 >      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
2207 >      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
2208 >      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
2209  
2210        counter++;      
2211      }
# Line 1824 | Line 2248 | void OSUAnalysis::fill1DHistogram(TH1* h
2248    bool sameObjects = false;
2249    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2250  
2251 <  int counter = 0;  
2251 >  int pairCounter = 0;  
2252    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2253      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2254        
# Line 1833 | Line 2257 | void OSUAnalysis::fill1DHistogram(TH1* h
2257        //only take objects which have passed all cuts and pairs which have passed all cuts
2258        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2259        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2260 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue;
2260 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2261  
2262        string currentString = parameters.inputVariables.at(0);
2263        string inputVariable = "";
# Line 1850 | Line 2274 | void OSUAnalysis::fill1DHistogram(TH1* h
2274        double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2275        histo->Fill(value,puScaleFactor);
2276  
2277 <      counter++;      
2277 >      pairCounter++;      
2278      }
2279    }
2280  
# Line 1903 | Line 2327 | void OSUAnalysis::fill2DHistogram(TH2* h
2327    bool sameObjects = false;
2328    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2329  
2330 <  int counter = 0;  
2330 >  int pairCounter = 0;  
2331    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2332      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2333        
# Line 1912 | Line 2336 | void OSUAnalysis::fill2DHistogram(TH2* h
2336        //only take objects which have passed all cuts and pairs which have passed all cuts
2337        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2338        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2339 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue;
2339 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2340  
2341        string currentString = parameters.inputVariables.at(0);
2342        string inputVariable = "";
# Line 1943 | Line 2367 | void OSUAnalysis::fill2DHistogram(TH2* h
2367  
2368        histo->Fill(valueX,valueY,puScaleFactor);
2369  
2370 <      counter++;      
2370 >      pairCounter++;      
2371      }
2372    }
2373  
2374   }
2375  
2376  
2377 + template <class InputObject>
2378 + int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
2379 +
2380 +  int bestMatchIndex = -1;
2381 +  double bestMatchDeltaR = 999;
2382 +
2383 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2384 +
2385 +    double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
2386 +    if(currentDeltaR > 0.05) continue;
2387 +
2388 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2389 +      bestMatchIndex = mcparticle - mcparticles->begin();
2390 +      bestMatchDeltaR = currentDeltaR;
2391 +    }
2392 +
2393 +  }
2394 +
2395 +  return bestMatchIndex;
2396 +
2397 + }
2398 +
2399 +
2400 + int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
2401 +
2402 +  int bestMatchIndex = -1;
2403 +  double bestMatchDeltaR = 999;
2404 +
2405 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2406 +
2407 +    if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
2408 +
2409 +    double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
2410 +    if(currentDeltaR > 0.05) continue;
2411 +
2412 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2413 +      bestMatchIndex = mcparticle - mcparticles->begin();
2414 +      bestMatchDeltaR = currentDeltaR;
2415 +    }
2416 +
2417 +  }
2418 +
2419 +  return bestMatchIndex;
2420 + }
2421 +
2422 +
2423 + // bin      particle type
2424 + // ---      -------------
2425 + //  0        unmatched
2426 + //  1        u
2427 + //  2        d
2428 + //  3        s
2429 + //  4        c
2430 + //  5        b
2431 + //  6        t
2432 + //  7        e
2433 + //  8        mu
2434 + //  9        tau
2435 + // 10        nu
2436 + // 11        g
2437 + // 12        gamma
2438 + // 13        Z
2439 + // 14        W
2440 + // 15        light meson
2441 + // 16        K meson
2442 + // 17        D meson
2443 + // 18        B meson
2444 + // 19        light baryon
2445 + // 20        strange baryon
2446 + // 21        charm baryon
2447 + // 22        bottom baryon
2448 + // 23        other
2449 +
2450 + int OSUAnalysis::getPdgIdBinValue(int pdgId){
2451 +
2452 +  int binValue = -999;
2453 +  int absPdgId = fabs(pdgId);  
2454 +  if(pdgId == -1) binValue = 0;
2455 +  else if(absPdgId <= 6 ) binValue = absPdgId;
2456 +  else if(absPdgId == 11 ) binValue = 7;
2457 +  else if(absPdgId == 13 ) binValue = 8;
2458 +  else if(absPdgId == 15 ) binValue = 9;
2459 +  else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
2460 +  else if(absPdgId == 21 ) binValue = 11;
2461 +  else if(absPdgId == 22 ) binValue = 12;
2462 +  else if(absPdgId == 23 ) binValue = 13;
2463 +  else if(absPdgId == 24 ) binValue = 14;
2464 +  else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130  ) binValue = 15;
2465 +  else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400)  ) binValue = 16;
2466 +  else if(absPdgId > 400 && absPdgId < 500  ) binValue = 17;
2467 +  else if(absPdgId > 500 && absPdgId < 600  ) binValue = 18;
2468 +  else if(absPdgId > 1000 && absPdgId < 3000  ) binValue = 19;
2469 +  else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
2470 +  else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
2471 +  else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
2472 +
2473 +  else binValue = 23;
2474 +
2475 +  return binValue;
2476 +
2477 + }
2478 +
2479  
2480   DEFINE_FWK_MODULE(OSUAnalysis);
2481  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines