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.24 by lantonel, Mon Mar 4 23:25:13 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  
44 +  //always get the MC particles to do GEN-matching
45 +  objectsToGet.push_back("mcparticles");
46 +
47    //always get the event collection to do pile-up reweighting
48 <  allNecessaryObjects.push_back("events");
48 >  objectsToGet.push_back("events");
49  
50    //parse the histogram definitions
51    for(uint currentHistogramSet = 0; currentHistogramSet != histogramSets_.size(); currentHistogramSet++){
# Line 50 | Line 53 | OSUAnalysis::OSUAnalysis (const edm::Par
53      string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection");
54      if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs";
55      if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
56 +      objectsToGet.push_back(tempInputCollection);
57        objectsToPlot.push_back(tempInputCollection);
58 <      allNecessaryObjects.push_back(tempInputCollection);
58 >      objectsToCut.push_back(tempInputCollection);
59      }
60 <    else{//pair of objects, need to add them both to the things to allNecessaryObjects
60 >    else{//pair of objects, need to add them both to the things to objectsToGet
61        int dashIndex = tempInputCollection.find("-");
62        int spaceIndex = tempInputCollection.find(" ");
63        int secondWordLength = spaceIndex - dashIndex;
64 <      allNecessaryObjects.push_back(tempInputCollection);
65 <      allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s");
66 <      allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
64 >      objectsToGet.push_back(tempInputCollection);
65 >      objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
66 >      objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
67        objectsToPlot.push_back(tempInputCollection);
68        objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s");
69        objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
70 +      objectsToCut.push_back(tempInputCollection);
71 +      objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
72 +      objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
73 +
74        }
75  
76      vector<edm::ParameterSet> histogramList_  (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms"));
# Line 80 | Line 88 | OSUAnalysis::OSUAnalysis (const edm::Par
88  
89      }
90    }
91 <  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object
91 >  //make unique vector of objects we need to plot (so we can book a histogram with the number of each object)
92    sort( objectsToPlot.begin(), objectsToPlot.end() );
93    objectsToPlot.erase( unique( objectsToPlot.begin(), objectsToPlot.end() ), objectsToPlot.end() );
94  
95  
96  
97 +  //add histograms with the gen-matched id, mother id, and grandmother id
98 +  for(uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
99 +
100 +    string currentObject = objectsToPlot.at(currentObjectIndex);
101 +    if(currentObject != "muons" && currentObject != "electrons" && currentObject != "taus" && currentObject != "tracks" && currentObject != "photons" && currentObject != "superclusters") continue;
102 +
103 +    histogram tempIdHisto;
104 +    histogram tempMomIdHisto;
105 +    histogram tempGmaIdHisto;
106 +
107 +    tempIdHisto.inputCollection = currentObject;
108 +    tempMomIdHisto.inputCollection = currentObject;
109 +    tempGmaIdHisto.inputCollection = currentObject;
110 +
111 +    currentObject = currentObject.substr(0, currentObject.size()-1);
112 +    tempIdHisto.name = currentObject+"GenMatchId";
113 +    tempMomIdHisto.name = currentObject+"GenMatchMotherId";
114 +    tempGmaIdHisto.name = currentObject+"GenMatchGrandmotherId";
115 +
116 +    currentObject.at(0) = toupper(currentObject.at(0));
117 +    tempIdHisto.title = currentObject+" Gen-matched Particle";
118 +    tempMomIdHisto.title = currentObject+" Gen-matched Particle's Mother";
119 +    tempGmaIdHisto.title = currentObject+" Gen-matched Particle's Grandmother";
120 +
121 +    int maxNum = 24;
122 +    vector<double> binVector;
123 +    binVector.push_back(maxNum);
124 +    binVector.push_back(0);
125 +    binVector.push_back(maxNum);
126 +
127 +    tempIdHisto.bins = binVector;
128 +    tempIdHisto.inputVariables.push_back("genMatchedId");
129 +    tempMomIdHisto.bins = binVector;
130 +    tempMomIdHisto.inputVariables.push_back("genMatchedMotherId");
131 +    tempGmaIdHisto.bins = binVector;
132 +    tempGmaIdHisto.inputVariables.push_back("genMatchedGrandmotherId");
133 +
134 +    histograms.push_back(tempIdHisto);
135 +    histograms.push_back(tempMomIdHisto);
136 +    histograms.push_back(tempGmaIdHisto);
137 +    
138 +  }
139 +
140 +
141  
142    channel tempChannel;
143    //loop over all channels (event selections)
# Line 104 | Line 156 | OSUAnalysis::OSUAnalysis (const edm::Par
156        triggerNames   = channels_.at(currentChannel).getParameter<vector<string> >("triggers");
157        for(uint trigger = 0; trigger!= triggerNames.size(); trigger++)
158          tempChannel.triggers.push_back(triggerNames.at(trigger));
159 <      allNecessaryObjects.push_back("triggers");
159 >      objectsToGet.push_back("triggers");
160      }
161  
162  
# Line 146 | Line 198 | OSUAnalysis::OSUAnalysis (const edm::Par
198        }
199  
200  
201 +      if(currentHistogram.name.find("GenMatch")==std::string::npos) continue;
202 +
203 + // bin      particle type
204 + // ---      -------------
205 + //  0        unmatched
206 + //  1        u
207 + //  2        d
208 + //  3        s
209 + //  4        c
210 + //  5        b
211 + //  6        t
212 + //  7        e
213 + //  8        mu
214 + //  9        tau
215 + // 10        nu
216 + // 11        g
217 + // 12        gamma
218 + // 13        Z
219 + // 14        W
220 + // 15        light meson
221 + // 16        K meson
222 + // 17        D meson
223 + // 18        B meson
224 + // 19        light baryon
225 + // 20        strange baryon
226 + // 21        charm baryon
227 + // 22        bottom baryon
228 + // 23        other
229 +
230 +      vector<TString> labelArray;
231 +      labelArray.push_back("unmatched");
232 +      labelArray.push_back("u");
233 +      labelArray.push_back("d");
234 +      labelArray.push_back("s");
235 +      labelArray.push_back("c");
236 +      labelArray.push_back("b");
237 +      labelArray.push_back("t");
238 +      labelArray.push_back("e");
239 +      labelArray.push_back("#mu");
240 +      labelArray.push_back("#tau");
241 +      labelArray.push_back("#nu");
242 +      labelArray.push_back("g");
243 +      labelArray.push_back("#gamma");
244 +      labelArray.push_back("Z");
245 +      labelArray.push_back("W");
246 +      labelArray.push_back("light meson");
247 +      labelArray.push_back("K meson");
248 +      labelArray.push_back("D meson");
249 +      labelArray.push_back("B meson");
250 +      labelArray.push_back("light baryon");
251 +      labelArray.push_back("strange baryon");
252 +      labelArray.push_back("charm baryon");
253 +      labelArray.push_back("bottom baryon");
254 +      labelArray.push_back("other");
255 +
256 +      for(int bin = 0; bin !=currentHistogram.bins.at(0); bin++){
257 +        oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin));
258 +      }
259      }
260      //book a histogram for the number of each object type to be plotted
261      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
# Line 171 | Line 281 | OSUAnalysis::OSUAnalysis (const edm::Par
281      }
282  
283  
284 +    
285 +
286 +
287      //get list of cuts for this channel
288      vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
289  
# Line 181 | Line 294 | OSUAnalysis::OSUAnalysis (const edm::Par
294        string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
295        tempCut.inputCollection = tempInputCollection;
296        if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object
297 <        allNecessaryObjects.push_back(tempInputCollection);
297 >        objectsToGet.push_back(tempInputCollection);
298 >        objectsToCut.push_back(tempInputCollection);
299        }
300 <      else{//pair of objects, need to add them both to the things to allNecessaryObjects
300 >      else{//pair of objects, need to add them both to the things to objectsToGet
301          int dashIndex = tempInputCollection.find("-");
302          int spaceIndex = tempInputCollection.find(" ");
303          int secondWordLength = spaceIndex - dashIndex;
304 <        allNecessaryObjects.push_back(tempInputCollection);
305 <        allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s");
306 <        allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
304 >        objectsToGet.push_back(tempInputCollection);
305 >        objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s");
306 >        objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
307 >        objectsToCut.push_back(tempInputCollection);
308 >        objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s");
309 >        objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s");
310  
311        }
312  
# Line 259 | Line 376 | OSUAnalysis::OSUAnalysis (const edm::Par
376    }//end loop over channels
377  
378  
379 <  //make unique vector of objects we need to cut on (so we make sure to get them from the event)
380 <  sort( allNecessaryObjects.begin(), allNecessaryObjects.end() );
381 <  allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() );
379 >  //make unique vector of objects we need to get from the event
380 >  sort( objectsToGet.begin(), objectsToGet.end() );
381 >  objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() );
382 >  //make unique vector of objects we need to cut on
383 >  sort( objectsToCut.begin(), objectsToCut.end() );
384 >  objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() );
385  
386  
387   }
# Line 279 | Line 399 | void
399   OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
400   {
401  
282
402    // Retrieve necessary collections from the event.
403 <  edm::Handle<BNtriggerCollection> triggers;
404 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "triggers") != allNecessaryObjects.end())
403 >
404 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end())
405      event.getByLabel (triggers_, triggers);
406  
407 <  edm::Handle<BNjetCollection> jets;
289 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end())
407 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end())
408      event.getByLabel (jets_, jets);
409  
410 <  edm::Handle<BNmuonCollection> muons;
293 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end())
410 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end())
411      event.getByLabel (muons_, muons);
412  
413 <  edm::Handle<BNelectronCollection> electrons;
297 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end())
413 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end())
414      event.getByLabel (electrons_, electrons);
415  
416 <  edm::Handle<BNeventCollection> events;
301 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end())
416 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end())
417      event.getByLabel (events_, events);
418  
419 <  edm::Handle<BNtauCollection> taus;
305 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end())
419 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end())
420      event.getByLabel (taus_, taus);
421  
422 <  edm::Handle<BNmetCollection> mets;
309 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end())
422 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end())
423      event.getByLabel (mets_, mets);
424  
425 <  edm::Handle<BNtrackCollection> tracks;
313 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end())
425 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end())
426      event.getByLabel (tracks_, tracks);
427  
428 <  edm::Handle<BNgenjetCollection> genjets;
317 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end())
428 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end())
429      event.getByLabel (genjets_, genjets);
430  
431 <  edm::Handle<BNmcparticleCollection> mcparticles;
321 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end())
431 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end())
432      event.getByLabel (mcparticles_, mcparticles);
433  
434 <  edm::Handle<BNprimaryvertexCollection> primaryvertexs;
325 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end())
434 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end())
435      event.getByLabel (primaryvertexs_, primaryvertexs);
436  
437 <  edm::Handle<BNbxlumiCollection> bxlumis;
329 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end())
437 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end())
438      event.getByLabel (bxlumis_, bxlumis);
439  
440 <  edm::Handle<BNphotonCollection> photons;
333 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end())
440 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end())
441      event.getByLabel (photons_, photons);
442  
443 <  edm::Handle<BNsuperclusterCollection> superclusters;
337 <  if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end())
443 >  if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end())
444      event.getByLabel (superclusters_, superclusters);
445  
446  
447    //get pile-up event weight
448 <  double puScaleFactor = 0;
448 >  double puScaleFactor = 1.00;
449    if(datasetType_ != "data")
450      puScaleFactor = puWeight_->at (events->at (0).numTruePV);
345  else
346    puScaleFactor = 1.00;
347
451  
452    //loop over all channels
453  
# Line 362 | Line 465 | OSUAnalysis::analyze (const edm::Event &
465      }
466  
467      //loop over all cuts
468 +
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 406 | Line 527 | OSUAnalysis::analyze (const edm::Event &
527      eventPassedAllCuts = eventPassedAllCuts && triggerDecision;
528  
529  
530 +
531      for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){
532  
533        //loop over all objects and count how many passed the cumulative selection up to this point
# Line 430 | Line 552 | OSUAnalysis::analyze (const edm::Event &
552      if(!eventPassedAllCuts)continue;
553  
554  
555 <
555 >    
556  
557      //set position of primary vertex in event, in order to calculate quantities relative to it
558 <    primaryVertex_ = 0;
559 <    vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
560 <    for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
561 <      if(!vertexFlags.at(vertexIndex)) continue;
562 <      primaryVertex_ = new BNprimaryvertex (primaryvertexs->at (vertexIndex));
563 <      break;
558 >    if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) {
559 >      vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back();
560 >      for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){
561 >        if(!vertexFlags.at(vertexIndex)) continue;
562 >        chosenPrimaryVertex = & primaryvertexs->at(vertexIndex);
563 >        break;
564 >      }
565 >    }
566 >    else {
567 >      chosenPrimaryVertex = & primaryvertexs->at(0);
568      }
569  
570  
# Line 450 | Line 576 | OSUAnalysis::analyze (const edm::Event &
576        if(currentHistogram.inputVariables.size() == 1){
577          TH1D* histo;
578          histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name);
579 +
580 +
581 +
582          if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
583 <        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
583 >        else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
584          else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),\
585                                                                                         cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(),\
586                                                                                         cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor);
# Line 476 | Line 605 | OSUAnalysis::analyze (const edm::Event &
605        else if(currentHistogram.inputVariables.size() == 2){
606          TH2D* histo;
607          histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name);
608 +
609 +
610 +
611          if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor);
612          else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor);
613          else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \
# Line 501 | Line 633 | OSUAnalysis::analyze (const edm::Event &
633        }
634      }
635  
504
636      //fills histograms with the sizes of collections
637      for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){
638        string currentObject = objectsToPlot.at(currentObjectIndex);
639  
640 <      if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs";
641 <      else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs";
642 <      else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs";
643 <      string tempCurrentObject = currentObject;
640 >      string objectToPlot = "";
641 >
642 >      if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs";
643 >      else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs";
644 >      else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs";
645 >      else objectToPlot = currentObject;
646 >      string tempCurrentObject = objectToPlot;
647        tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0));
648        string histoName = "num" + tempCurrentObject;
649  
650  
651 <      if(currentObject == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
652 <      else if(currentObject == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
653 <      else if(currentObject == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
654 <      else if(currentObject == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
655 <      else if(currentObject == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
656 <      else if(currentObject == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
657 <      else if(currentObject == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
658 <      else if(currentObject == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
659 <      else if(currentObject == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
660 <      else if(currentObject == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
661 <      else if(currentObject == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
662 <      else if(currentObject == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
663 <      else if(currentObject == "primaryvertexs"){
651 >
652 >
653 >      //set position of primary vertex in event, in order to calculate quantities relative to it
654 >      if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) {
655 >        vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back();
656 >        int numToPlot = 0;
657 >        for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){
658 >          if(lastCutFlags.at(currentFlag)) numToPlot++;
659 >        }
660 >        if(objectToPlot == "primaryvertexs"){
661 >          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
662 >          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
663 >        }
664 >        else
665 >          oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,puScaleFactor);        
666 >      }
667 >      else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor);
668 >      else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor);
669 >      else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor);
670 >      else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor);
671 >      else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor);
672 >      else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor);
673 >      else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor);
674 >      else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor);
675 >      else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor);
676 >      else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor);
677 >      else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor);
678 >      else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor);
679 >      else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor);
680 >      else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor);
681 >      else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor);
682 >      else if(objectToPlot == "primaryvertexs"){
683          oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size());
684          oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor);
685        }
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);
686  
687      }
688 +
689 +
690 +
691 +
692    } //end loop over channel
693  
694    masterCutFlow_->fillCutFlow(puScaleFactor);
# Line 736 | Line 890 | OSUAnalysis::valueLookup (const BNmuon*
890    else if(variable == "ecalIso") value = object->ecalIso;
891    else if(variable == "hcalIso") value = object->hcalIso;
892    else if(variable == "caloIso") value = object->caloIso;
893 <  else if(variable == "trackIsoDR03") value = object->trackIsoDR03;
893 >  else if(variable == "trackIsDR03") value = object->trackIsoDR03;
894    else if(variable == "ecalIsoDR03") value = object->ecalIsoDR03;
895    else if(variable == "hcalIsoDR03") value = object->hcalIsoDR03;
896    else if(variable == "caloIsoDR03") value = object->caloIsoDR03;
# Line 875 | Line 1029 | OSUAnalysis::valueLookup (const BNmuon*
1029    else if(variable == "time_ndof") value = object->time_ndof;
1030  
1031    //user-defined variables
1032 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1033 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1034 <  else if(variable == "detIso") value = (object->trackIso) / object->pt;
1032 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1033 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1034 >  else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt;
1035    else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt;
1036    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt;
1037    else if(variable == "tightID") {
# Line 888 | Line 1042 | OSUAnalysis::valueLookup (const BNmuon*
1042        && object->numberOfMatchedStations > 1    \
1043        && fabs(object->correctedD0Vertex) < 0.2  \
1044        && fabs(object->correctedDZ) < 0.5        \
1045 <      && object->numberOfValidPixelHits > 0     \
1045 >      && object->numberOfValidPixelHits > 0             \
1046        && object->numberOfLayersWithMeasurement > 5;
893    
894
1047    }
1048    else if(variable == "tightIDdisplaced"){
1049      value = object->isGlobalMuon > 0            \
898      && object->isPFMuon > 0                   \
1050        && object->normalizedChi2 < 10            \
1051        && object->numberOfValidMuonHits > 0      \
1052        && object->numberOfMatchedStations > 1    \
# Line 903 | Line 1054 | OSUAnalysis::valueLookup (const BNmuon*
1054        && object->numberOfLayersWithMeasurement > 5;
1055    }
1056  
1057 +  else if(variable == "genMatchedId"){
1058 +    int index = getGenMatchedParticleIndex(object);
1059 +    if(index == -1) value = 0;
1060 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1061 +  }
1062 +  else if(variable == "genMatchedMotherId"){
1063 +    int index = getGenMatchedParticleIndex(object);
1064 +    if(index == -1) value = 0;
1065 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1066 +  }
1067 +  else if(variable == "genMatchedGrandmotherId"){
1068 +    int index = getGenMatchedParticleIndex(object);
1069 +    if(index == -1) value = 0;
1070 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1071 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1072 +      if(motherIndex == -1) value = 0;
1073 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1074 +    }
1075 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1076 +  }
1077 +
1078 +
1079 +
1080    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1081  
1082    value = applyFunction(function, value);
# Line 1073 | Line 1247 | OSUAnalysis::valueLookup (const BNelectr
1247    else if(variable == "passConvVeto") value = object->passConvVeto;
1248  
1249    //user-defined variables
1250 <  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1251 <  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (primaryVertex_->xError, primaryVertex_->yError));
1250 >  else if(variable == "correctedD0VertexErr") value =  hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1251 >  else if(variable == "correctedD0VertexSig") value =  object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError));
1252    else if(variable == "detIso") value = (object->trackIso) / object->pt;
1253    else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt;
1254 +  else if(variable == "correctedD0VertexInEB"){
1255 +    if(fabs(object->eta) < 0.8) value = object->correctedD0Vertex;
1256 +    else value = -999;
1257 +  }
1258 +  else if(variable == "correctedD0VertexOutEB"){
1259 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0Vertex;
1260 +    else value = -999;
1261 +  }
1262 +  else if(variable == "correctedD0VertexEE"){
1263 +    if(object->isEE) value = object->correctedD0Vertex;
1264 +    else value = -999;
1265 +  }
1266 +
1267 +  else if(variable == "correctedD0BeamspotInEB"){
1268 +    if(fabs(object->eta) < 0.8) value = object->correctedD0;
1269 +    else value = -999;
1270 +  }
1271 +  else if(variable == "correctedD0BeamspotOutEB"){
1272 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0;
1273 +    else value = -999;
1274 +  }
1275 +  else if(variable == "correctedD0BeamspotEE"){
1276 +    if(object->isEE) value = object->correctedD0;
1277 +    else value = -999;
1278 +  }
1279 +
1280 +  else if(variable == "correctedD0OriginInEB"){
1281 +    if(fabs(object->eta) < 0.8) value = object->tkD0;
1282 +    else value = -999;
1283 +  }
1284 +  else if(variable == "correctedD0OriginOutEB"){
1285 +    if(object->isEB && fabs(object->eta) > 0.8) value = object->tkD0;
1286 +    else value = -999;
1287 +  }
1288 +  else if(variable == "correctedD0OriginEE"){
1289 +    if(object->isEE) value = object->tkD0;
1290 +    else value = -999;
1291 +  }
1292 +
1293 +  else if(variable == "tightIDdisplaced"){
1294 +    if (object->isEB)
1295 +      {
1296 +        value = fabs(object->delEtaIn) < 0.004 \
1297 +          && fabs (object->delPhiIn) < 0.03 \
1298 +          && object->scSigmaIEtaIEta < 0.01 \
1299 +          && object->hadOverEm < 0.12 \
1300 +          && object->absInvEMinusInvPin < 0.05;
1301 +      }
1302 +    else
1303 +      {
1304 +        value = fabs (object->delEtaIn) < 0.005 \
1305 +          && fabs (object->delPhiIn) < 0.02 \
1306 +          && object->scSigmaIEtaIEta < 0.03 \
1307 +          && object->hadOverEm < 0.10 \
1308 +          && object->absInvEMinusInvPin < 0.05;
1309 +      }
1310 +  }
1311 +
1312 +  else if(variable == "genMatchedId"){
1313 +    int index = getGenMatchedParticleIndex(object);
1314 +    if(index == -1) value = 0;
1315 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1316 +  }
1317 +  else if(variable == "genMatchedMotherId"){
1318 +    int index = getGenMatchedParticleIndex(object);
1319 +    if(index == -1) value = 0;
1320 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1321 +  }
1322 +  else if(variable == "genMatchedGrandmotherId"){
1323 +    int index = getGenMatchedParticleIndex(object);
1324 +    if(index == -1) value = 0;
1325 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1326 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1327 +      if(motherIndex == -1) value = 0;
1328 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1329 +    }
1330 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1331 +  }
1332 +
1333  
1334  
1335    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
# Line 1210 | Line 1463 | OSUAnalysis::valueLookup (const BNtau* o
1463    else if(variable == "leadingTrackValid") value = object->leadingTrackValid;
1464  
1465  
1466 +  else if(variable == "genMatchedId"){
1467 +    int index = getGenMatchedParticleIndex(object);
1468 +    if(index == -1) value = 0;
1469 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1470 +  }
1471 +  else if(variable == "genMatchedMotherId"){
1472 +    int index = getGenMatchedParticleIndex(object);
1473 +    if(index == -1) value = 0;
1474 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1475 +  }
1476 +  else if(variable == "genMatchedGrandmotherId"){
1477 +    int index = getGenMatchedParticleIndex(object);
1478 +    if(index == -1) value = 0;
1479 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1480 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1481 +      if(motherIndex == -1) value = 0;
1482 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1483 +    }
1484 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1485 +  }
1486 +
1487 +
1488    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1489  
1490    value = applyFunction(function, value);
# Line 1313 | Line 1588 | OSUAnalysis::valueLookup (const BNtrack*
1588    else if(variable == "isHighPurity") value = object->isHighPurity;
1589  
1590  
1591 +  else if(variable == "genMatchedId"){
1592 +    int index = getGenMatchedParticleIndex(object);
1593 +    if(index == -1) value = 0;
1594 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1595 +  }
1596 +  else if(variable == "genMatchedMotherId"){
1597 +    int index = getGenMatchedParticleIndex(object);
1598 +    if(index == -1) value = 0;
1599 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1600 +  }
1601 +  else if(variable == "genMatchedGrandmotherId"){
1602 +    int index = getGenMatchedParticleIndex(object);
1603 +    if(index == -1) value = 0;
1604 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1605 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1606 +      if(motherIndex == -1) value = 0;
1607 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1608 +    }
1609 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1610 +  }
1611 +
1612 +
1613 +
1614    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1615  
1616    value = applyFunction(function, value);
# Line 1440 | Line 1738 | OSUAnalysis::valueLookup (const BNmcpart
1738  
1739    //user-defined variables
1740    else if (variable == "d0"){
1741 <    double vx = object->vx - primaryVertex_->x,
1742 <      vy = object->vy - primaryVertex_->y,
1741 >    double vx = object->vx - chosenPrimaryVertex->x,
1742 >      vy = object->vy - chosenPrimaryVertex->y,
1743        px = object->px,
1744        py = object->py,
1745        pt = object->pt;
1746      value = (-vx * py + vy * px) / pt;
1747    }
1748    else if (variable == "dz"){
1749 <    double vx = object->vx - primaryVertex_->x,
1750 <      vy = object->vy - primaryVertex_->y,
1751 <      vz = object->vz - primaryVertex_->z,
1749 >    double vx = object->vx - chosenPrimaryVertex->x,
1750 >      vy = object->vy - chosenPrimaryVertex->y,
1751 >      vz = object->vz - chosenPrimaryVertex->z,
1752        px = object->px,
1753        py = object->py,
1754        pz = object->pz,
1755        pt = object->pt;
1756      value = vz - (vx * px + vy * py)/pt * (pz/pt);
1757    }
1758 <
1758 >  else if(variable == "v0"){
1759 >    value = sqrt(object->vx*object->vx + object->vy*object->vy);
1760 >  }
1761 >  else if(variable == "deltaV0"){
1762 >    value = sqrt((object->vx-chosenPrimaryVertex->x)*(object->vx-chosenPrimaryVertex->x) + (object->vy-chosenPrimaryVertex->y)*(object->vy-chosenPrimaryVertex->y));
1763 >  }
1764 >  else if (variable == "deltaVx"){
1765 >    value = object->vx - chosenPrimaryVertex->x;
1766 >  }
1767 >  else if (variable == "deltaVy"){
1768 >    value = object->vy - chosenPrimaryVertex->y;
1769 >  }
1770 >  else if (variable == "deltaVz"){
1771 >    value = object->vz - chosenPrimaryVertex->z;
1772 >  }
1773  
1774  
1775    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
# Line 1587 | Line 1899 | OSUAnalysis::valueLookup (const BNphoton
1899    else if(variable == "seedRecoFlag") value = object->seedRecoFlag;
1900  
1901  
1902 +  else if(variable == "genMatchedId"){
1903 +    int index = getGenMatchedParticleIndex(object);
1904 +    if(index == -1) value = 0;
1905 +    else value = getPdgIdBinValue(mcparticles->at(index).id);
1906 +  }
1907 +  else if(variable == "genMatchedMotherId"){
1908 +    int index = getGenMatchedParticleIndex(object);
1909 +    if(index == -1) value = 0;
1910 +    else value = getPdgIdBinValue(mcparticles->at(index).motherId);
1911 +  }
1912 +  else if(variable == "genMatchedGrandmotherId"){
1913 +    int index = getGenMatchedParticleIndex(object);
1914 +    if(index == -1) value = 0;
1915 +    else if(fabs(mcparticles->at(index).motherId) == 15){
1916 +      int motherIndex = findTauMotherIndex(&mcparticles->at(index));
1917 +      if(motherIndex == -1) value = 0;
1918 +      else value = getPdgIdBinValue(mcparticles->at(motherIndex).motherId);
1919 +    }
1920 +    else value = getPdgIdBinValue(mcparticles->at(index).grandMotherId);
1921 +  }
1922 +
1923 +
1924    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1925  
1926    value = applyFunction(function, value);
# Line 1622 | Line 1956 | OSUAnalysis::valueLookup (const BNmuon*
1956  
1957    double value = 0.0;
1958  
1959 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
1959 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
1960 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
1961    else if(variable == "invMass"){
1962      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
1963      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
1964 <    value = (fourVector1 + fourVector2).M();
1965 <  }
1964 >
1965 >    value = (fourVector1 + fourVector2).M();}
1966 >  else if(variable == "threeDAngle")
1967 >    {
1968 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
1969 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
1970 >      value = (threeVector1.Angle(threeVector2));
1971 >    }
1972 >
1973 >
1974 >
1975    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
1976    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
1977 <  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
1978 <
1977 >  else if(variable == "d0Sign"){
1978 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
1979 >    if(d0Sign < 0) value = -0.5;
1980 >    else if (d0Sign > 0) value = 0.5;
1981 >    else value = -999;
1982 >  }
1983 >  else if(variable == "muon1CorrectedD0Vertex"){
1984 >    value = object1->correctedD0Vertex;
1985 >  }
1986 >  else if(variable == "muon2CorrectedD0Vertex"){
1987 >    value = object2->correctedD0Vertex;
1988 >  }
1989    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
1990  
1991 +
1992 +
1993 +
1994 +
1995 +
1996 +
1997 +
1998    value = applyFunction(function, value);
1999  
2000    return value;
# Line 1644 | Line 2005 | OSUAnalysis::valueLookup (const BNelectr
2005  
2006    double value = 0.0;
2007  
2008 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
2008 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2009 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2010    else if(variable == "invMass"){
2011      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2012      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2013 <    value = (fourVector1 + fourVector2).M();
2014 <  }
2013 >
2014 >    value = (fourVector1 + fourVector2).M();}
2015 > else if(variable == "threeDAngle")
2016 >    {
2017 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2018 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2019 >      value = (threeVector1.Angle(threeVector2));
2020 >    }
2021 >  
2022 >
2023 >
2024 >
2025    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2026    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2027 +
2028    else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2029  
2030 +
2031 +
2032 +  else if(variable == "d0Sign"){
2033 +    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2034 +    if(d0Sign < 0) value = -0.5;
2035 +    else if (d0Sign > 0) value = 0.5;
2036 +    else value = -999;
2037 +  }
2038 +  else if(variable == "electron1CorrectedD0Vertex"){
2039 +    value = object1->correctedD0Vertex;
2040 +  }
2041 +  else if(variable == "electron2CorrectedD0Vertex"){
2042 +    value = object2->correctedD0Vertex;
2043 +  }
2044 +
2045    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2046  
2047    value = applyFunction(function, value);
# Line 1666 | Line 2054 | OSUAnalysis::valueLookup (const BNelectr
2054  
2055    double value = 0.0;
2056  
2057 <  if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi);
2057 >  if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi));
2058 >  else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi);
2059    else if(variable == "invMass"){
2060      TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy);
2061      TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy);
2062 <    value = (fourVector1 + fourVector2).M();
2063 <  }
2062 >
2063 >    value = (fourVector1 + fourVector2).M();}
2064 > else if(variable == "threeDAngle")
2065 >    {
2066 >      TVector3 threeVector1(object1->px, object1->py, object1->pz);
2067 >      TVector3 threeVector2(object2->px, object2->py, object2->pz);
2068 >      value = (threeVector1.Angle(threeVector2));
2069 >    }
2070 >
2071 >
2072    else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex;
2073    else if(variable == "deltaAbsCorrectedD0Vertex") value = fabs(object1->correctedD0Vertex) - fabs(object2->correctedD0Vertex);
2074 <  else if(variable == "d0Sign") value = object1->correctedD0Vertex*object2->correctedD0Vertex/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2074 >  else if(variable == "d0Sign"){
2075 >    double d0Sign = (object1->correctedD0Vertex*object2->correctedD0Vertex)/fabs(object1->correctedD0Vertex*object2->correctedD0Vertex);
2076 >    if(d0Sign < 0) value = -0.5;
2077 >    else if (d0Sign > 0) value = 0.5;
2078 >    else value = -999;
2079 >  }
2080 >  else if(variable == "electronCorrectedD0Vertex"){
2081 >    value = object1->correctedD0Vertex;
2082 >  }
2083 >  else if(variable == "muonCorrectedD0Vertex"){
2084 >    value = object2->correctedD0Vertex;
2085 >  }
2086 >  else if(variable == "electronCorrectedD0"){
2087 >    value = object1->correctedD0;
2088 >  }
2089 >  else if(variable == "muonCorrectedD0"){
2090 >    value = object2->correctedD0;
2091 >  }
2092 >
2093  
2094    else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
2095  
# Line 1688 | Line 2103 | double
2103   OSUAnalysis::applyFunction(string function, double value){
2104  
2105    if(function == "abs") value = fabs(value);
2106 +  else if(function == "fabs") value = fabs(value);
2107  
2108  
2109 +  else if(function == "") value = value;
2110 +  else{std::cout << "WARNING: invalid function '" << function << "'\n";}
2111 +
2112    return value;
2113  
2114   }
# Line 1698 | Line 2117 | OSUAnalysis::applyFunction(string functi
2117   template <class InputCollection>
2118   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){
2119  
2120 +
2121    for (uint object = 0; object != inputCollection->size(); object++){
2122  
2123 +
2124      bool decision = true;//object passes if this cut doesn't cut on that type of object
2125  
2126 +
2127      if(currentCut.inputCollection == inputType){
2128  
2129        vector<bool> subcutDecisions;
# Line 1723 | Line 2145 | void OSUAnalysis::setObjectFlags(cut &cu
2145      }
2146      individualFlags.at(inputType).at(currentCutIndex).push_back(decision);
2147  
2148 +
2149      //set flags for objects that pass each cut AND all the previous cuts
2150      bool previousCumulativeFlag = true;
2151      for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){
# Line 1731 | Line 2154 | void OSUAnalysis::setObjectFlags(cut &cu
2154      }
2155      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2156  
2157 +
2158    }
2159  
2160   }
# Line 1738 | Line 2162 | void OSUAnalysis::setObjectFlags(cut &cu
2162  
2163   template <class InputCollection1, class InputCollection2>
2164   void OSUAnalysis::setObjectFlags(cut &currentCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \
2165 <                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, string inputType){
2165 >                                 InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){
2166  
2167  
2168    bool sameObjects = false;
2169    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2170  
2171 +
2172    int counter = 0;  
2173    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2174      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
# Line 1781 | Line 2206 | void OSUAnalysis::setObjectFlags(cut &cu
2206          if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true;
2207          else{ previousCumulativeFlag = false; break;}
2208        }
2209 <      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision);
2209 >      //apply flags for the components of the composite object as well
2210 >      bool currentCumulativeFlag = true;
2211 >      if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision;
2212 >      else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2);
2213 >      else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1);
2214 >      else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2);
2215 >      cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag);
2216  
2217        counter++;      
2218      }
# Line 1824 | Line 2255 | void OSUAnalysis::fill1DHistogram(TH1* h
2255    bool sameObjects = false;
2256    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2257  
2258 <  int counter = 0;  
2258 >  int pairCounter = 0;  
2259    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2260      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2261        
# Line 1833 | Line 2264 | void OSUAnalysis::fill1DHistogram(TH1* h
2264        //only take objects which have passed all cuts and pairs which have passed all cuts
2265        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2266        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2267 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue;
2267 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2268  
2269        string currentString = parameters.inputVariables.at(0);
2270        string inputVariable = "";
# Line 1850 | Line 2281 | void OSUAnalysis::fill1DHistogram(TH1* h
2281        double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function);
2282        histo->Fill(value,puScaleFactor);
2283  
2284 <      counter++;      
2284 >      pairCounter++;      
2285      }
2286    }
2287  
# Line 1903 | Line 2334 | void OSUAnalysis::fill2DHistogram(TH2* h
2334    bool sameObjects = false;
2335    if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true;
2336  
2337 <  int counter = 0;  
2337 >  int pairCounter = 0;  
2338    for (uint object1 = 0; object1 != inputCollection1->size(); object1++){
2339      for (uint object2 = 0; object2 != inputCollection2->size(); object2++){
2340        
# Line 1912 | Line 2343 | void OSUAnalysis::fill2DHistogram(TH2* h
2343        //only take objects which have passed all cuts and pairs which have passed all cuts
2344        if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue;
2345        if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue;
2346 <      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue;
2346 >      if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue;
2347  
2348        string currentString = parameters.inputVariables.at(0);
2349        string inputVariable = "";
# Line 1943 | Line 2374 | void OSUAnalysis::fill2DHistogram(TH2* h
2374  
2375        histo->Fill(valueX,valueY,puScaleFactor);
2376  
2377 <      counter++;      
2377 >      pairCounter++;      
2378 >    }
2379 >  }
2380 >
2381 > }
2382 >
2383 >
2384 > template <class InputObject>
2385 > int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){
2386 >
2387 >  int bestMatchIndex = -1;
2388 >  double bestMatchDeltaR = 999;
2389 >
2390 >  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2391 >
2392 >    double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi);
2393 >    if(currentDeltaR > 0.05) continue;
2394 >
2395 >    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2396 >      bestMatchIndex = mcparticle - mcparticles->begin();
2397 >      bestMatchDeltaR = currentDeltaR;
2398      }
2399 +
2400    }
2401  
2402 +  return bestMatchIndex;
2403 +
2404   }
2405  
2406  
2407 + int OSUAnalysis::findTauMotherIndex(const BNmcparticle* tau){
2408 +
2409 +  int bestMatchIndex = -1;
2410 +  double bestMatchDeltaR = 999;
2411 +
2412 +  for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){
2413 +
2414 +    if(fabs(mcparticle->id) != 15 || mcparticle->status !=3) continue;
2415 +
2416 +    double currentDeltaR = deltaR(tau->eta,tau->phi,mcparticle->eta,mcparticle->phi);
2417 +    if(currentDeltaR > 0.05) continue;
2418 +
2419 +    if(currentDeltaR < bestMatchDeltaR && mcparticles->at(mcparticle - mcparticles->begin()).id != mcparticles->at(mcparticle - mcparticles->begin()).motherId){
2420 +      bestMatchIndex = mcparticle - mcparticles->begin();
2421 +      bestMatchDeltaR = currentDeltaR;
2422 +    }
2423 +
2424 +  }
2425 +
2426 +  return bestMatchIndex;
2427 + }
2428 +
2429 +
2430 + // bin      particle type
2431 + // ---      -------------
2432 + //  0        unmatched
2433 + //  1        u
2434 + //  2        d
2435 + //  3        s
2436 + //  4        c
2437 + //  5        b
2438 + //  6        t
2439 + //  7        e
2440 + //  8        mu
2441 + //  9        tau
2442 + // 10        nu
2443 + // 11        g
2444 + // 12        gamma
2445 + // 13        Z
2446 + // 14        W
2447 + // 15        light meson
2448 + // 16        K meson
2449 + // 17        D meson
2450 + // 18        B meson
2451 + // 19        light baryon
2452 + // 20        strange baryon
2453 + // 21        charm baryon
2454 + // 22        bottom baryon
2455 + // 23        other
2456 +
2457 + int OSUAnalysis::getPdgIdBinValue(int pdgId){
2458 +
2459 +  int binValue = -999;
2460 +  int absPdgId = fabs(pdgId);  
2461 +  if(pdgId == -1) binValue = 0;
2462 +  else if(absPdgId <= 6 ) binValue = absPdgId;
2463 +  else if(absPdgId == 11 ) binValue = 7;
2464 +  else if(absPdgId == 13 ) binValue = 8;
2465 +  else if(absPdgId == 15 ) binValue = 9;
2466 +  else if(absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) binValue = 10;
2467 +  else if(absPdgId == 21 ) binValue = 11;
2468 +  else if(absPdgId == 22 ) binValue = 12;
2469 +  else if(absPdgId == 23 ) binValue = 13;
2470 +  else if(absPdgId == 24 ) binValue = 14;
2471 +  else if(absPdgId > 100 && absPdgId < 300 && absPdgId != 130  ) binValue = 15;
2472 +  else if( absPdgId == 130 || (absPdgId > 300 && absPdgId < 400)  ) binValue = 16;
2473 +  else if(absPdgId > 400 && absPdgId < 500  ) binValue = 17;
2474 +  else if(absPdgId > 500 && absPdgId < 600  ) binValue = 18;
2475 +  else if(absPdgId > 1000 && absPdgId < 3000  ) binValue = 19;
2476 +  else if(absPdgId > 3000 && absPdgId < 4000  ) binValue = 20;
2477 +  else if(absPdgId > 4000 && absPdgId < 5000  ) binValue = 21;
2478 +  else if(absPdgId > 5000 && absPdgId < 6000  ) binValue = 22;
2479 +
2480 +  else binValue = 23;
2481 +
2482 +  return binValue;
2483 +
2484 + }
2485 +
2486  
2487   DEFINE_FWK_MODULE(OSUAnalysis);
2488  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines