17 |
|
superclusters_ (cfg.getParameter<edm::InputTag> ("superclusters")), |
18 |
|
triggers_ (cfg.getParameter<edm::InputTag> ("triggers")), |
19 |
|
puFile_ (cfg.getParameter<std::string> ("puFile")), |
20 |
+ |
deadEcalFile_ (cfg.getParameter<std::string> ("deadEcalFile")), |
21 |
|
dataPU_ (cfg.getParameter<std::string> ("dataPU")), |
22 |
|
dataset_ (cfg.getParameter<std::string> ("dataset")), |
23 |
|
datasetType_ (cfg.getParameter<std::string> ("datasetType")), |
30 |
|
TH1::SetDefaultSumw2 (); |
31 |
|
|
32 |
|
//create pile-up reweighting object, if necessary |
33 |
< |
if(datasetType_ != "data") puWeight_ = new PUWeight (puFile_, dataPU_, dataset_); |
33 |
> |
// if(datasetType_ != "data") puWeight_ = new PUWeight (puFile_, dataPU_, dataset_); |
34 |
|
|
35 |
|
// Construct Cutflow Objects. These store the results of cut decisions and |
36 |
|
// handle filling cut flow histograms. |
38 |
|
std::vector<TFileDirectory> directories; |
39 |
|
|
40 |
|
//always get vertex collection so we can assign the primary vertex in the event |
41 |
< |
allNecessaryObjects.push_back("primaryvertexs"); |
41 |
> |
|
42 |
> |
objectsToGet.push_back("primaryvertexs"); |
43 |
|
//always make the plot of number of primary verticex (to check pile-up reweighting) |
44 |
|
objectsToPlot.push_back("primaryvertexs"); |
45 |
|
|
46 |
|
//always get the MC particles to do GEN-matching |
47 |
< |
allNecessaryObjects.push_back("mcparticles"); |
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++){ |
55 |
|
string tempInputCollection = histogramSets_.at(currentHistogramSet).getParameter<string> ("inputCollection"); |
56 |
|
if(tempInputCollection == "muon-electron pairs") tempInputCollection = "electron-muon pairs"; |
57 |
|
if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object |
58 |
+ |
objectsToGet.push_back(tempInputCollection); |
59 |
|
objectsToPlot.push_back(tempInputCollection); |
60 |
< |
allNecessaryObjects.push_back(tempInputCollection); |
60 |
> |
objectsToCut.push_back(tempInputCollection); |
61 |
|
} |
62 |
< |
else{//pair of objects, need to add them both to the things to allNecessaryObjects |
62 |
> |
else{//pair of objects, need to add them both to the things to objectsToGet |
63 |
|
int dashIndex = tempInputCollection.find("-"); |
64 |
|
int spaceIndex = tempInputCollection.find(" "); |
65 |
|
int secondWordLength = spaceIndex - dashIndex; |
66 |
< |
allNecessaryObjects.push_back(tempInputCollection); |
67 |
< |
allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
68 |
< |
allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
66 |
> |
objectsToGet.push_back(tempInputCollection); |
67 |
> |
objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
68 |
> |
objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
69 |
|
objectsToPlot.push_back(tempInputCollection); |
70 |
|
objectsToPlot.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
71 |
|
objectsToPlot.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
72 |
+ |
objectsToCut.push_back(tempInputCollection); |
73 |
+ |
objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
74 |
+ |
objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
75 |
+ |
|
76 |
|
} |
77 |
|
|
78 |
|
vector<edm::ParameterSet> histogramList_ (histogramSets_.at(currentHistogramSet).getParameter<vector<edm::ParameterSet> >("histograms")); |
158 |
|
triggerNames = channels_.at(currentChannel).getParameter<vector<string> >("triggers"); |
159 |
|
for(uint trigger = 0; trigger!= triggerNames.size(); trigger++) |
160 |
|
tempChannel.triggers.push_back(triggerNames.at(trigger)); |
161 |
< |
allNecessaryObjects.push_back("triggers"); |
161 |
> |
objectsToGet.push_back("triggers"); |
162 |
|
} |
163 |
|
|
164 |
|
|
259 |
|
oneDHists_.at(currentChannel)[currentHistogram.name]->GetXaxis()->SetBinLabel(bin+1,labelArray.at(bin)); |
260 |
|
} |
261 |
|
} |
262 |
+ |
|
263 |
|
//book a histogram for the number of each object type to be plotted |
264 |
+ |
|
265 |
|
for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){ |
266 |
|
string currentObject = objectsToPlot.at(currentObjectIndex); |
267 |
|
int maxNum = 10; |
268 |
|
if(currentObject == "mcparticles") maxNum = 50; |
269 |
|
else if(currentObject == "primaryvertexs") maxNum = 50; |
270 |
< |
else if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs"; |
270 |
> |
else if(currentObject == "muon-muon pairs") |
271 |
> |
{ |
272 |
> |
currentObject = "dimuonPairs"; |
273 |
> |
} |
274 |
|
else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs"; |
275 |
|
else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs"; |
276 |
|
|
288 |
|
} |
289 |
|
|
290 |
|
|
279 |
– |
|
280 |
– |
|
291 |
|
|
292 |
|
|
293 |
|
|
301 |
|
string tempInputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection"); |
302 |
|
tempCut.inputCollection = tempInputCollection; |
303 |
|
if(tempInputCollection.find("pairs")==std::string::npos){ //just a single object |
304 |
< |
allNecessaryObjects.push_back(tempInputCollection); |
304 |
> |
objectsToGet.push_back(tempInputCollection); |
305 |
> |
objectsToCut.push_back(tempInputCollection); |
306 |
|
} |
307 |
< |
else{//pair of objects, need to add them both to the things to allNecessaryObjects |
307 |
> |
else{//pair of objects, need to add them both to the things to objectsToGet |
308 |
|
int dashIndex = tempInputCollection.find("-"); |
309 |
|
int spaceIndex = tempInputCollection.find(" "); |
310 |
|
int secondWordLength = spaceIndex - dashIndex; |
311 |
< |
allNecessaryObjects.push_back(tempInputCollection); |
312 |
< |
allNecessaryObjects.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
313 |
< |
allNecessaryObjects.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
311 |
> |
objectsToGet.push_back(tempInputCollection); |
312 |
> |
objectsToGet.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
313 |
> |
objectsToGet.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
314 |
> |
objectsToCut.push_back(tempInputCollection); |
315 |
> |
objectsToCut.push_back(tempInputCollection.substr(0,dashIndex)+"s"); |
316 |
> |
objectsToCut.push_back(tempInputCollection.substr(dashIndex+1,secondWordLength-1)+"s"); |
317 |
|
|
318 |
|
} |
319 |
|
|
383 |
|
}//end loop over channels |
384 |
|
|
385 |
|
|
386 |
< |
//make unique vector of objects we need to cut on (so we make sure to get them from the event) |
387 |
< |
sort( allNecessaryObjects.begin(), allNecessaryObjects.end() ); |
388 |
< |
allNecessaryObjects.erase( unique( allNecessaryObjects.begin(), allNecessaryObjects.end() ), allNecessaryObjects.end() ); |
386 |
> |
//make unique vector of objects we need to get from the event |
387 |
> |
sort( objectsToGet.begin(), objectsToGet.end() ); |
388 |
> |
objectsToGet.erase( unique( objectsToGet.begin(), objectsToGet.end() ), objectsToGet.end() ); |
389 |
> |
//make unique vector of objects we need to cut on |
390 |
> |
sort( objectsToCut.begin(), objectsToCut.end() ); |
391 |
> |
objectsToCut.erase( unique( objectsToCut.begin(), objectsToCut.end() ), objectsToCut.end() ); |
392 |
|
|
393 |
|
|
394 |
|
} |
406 |
|
OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup) |
407 |
|
{ |
408 |
|
|
392 |
– |
|
409 |
|
// Retrieve necessary collections from the event. |
410 |
|
|
411 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "triggers") != allNecessaryObjects.end()) |
411 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "triggers") != objectsToGet.end()) |
412 |
|
event.getByLabel (triggers_, triggers); |
413 |
|
|
414 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "jets") != allNecessaryObjects.end()) |
414 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "jets") != objectsToGet.end()) |
415 |
|
event.getByLabel (jets_, jets); |
416 |
|
|
417 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "muons") != allNecessaryObjects.end()) |
417 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "muons") != objectsToGet.end()) |
418 |
|
event.getByLabel (muons_, muons); |
419 |
|
|
420 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "electrons") != allNecessaryObjects.end()) |
420 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "electrons") != objectsToGet.end()) |
421 |
|
event.getByLabel (electrons_, electrons); |
422 |
|
|
423 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "events") != allNecessaryObjects.end()) |
423 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "events") != objectsToGet.end()) |
424 |
|
event.getByLabel (events_, events); |
425 |
|
|
426 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "taus") != allNecessaryObjects.end()) |
426 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "taus") != objectsToGet.end()) |
427 |
|
event.getByLabel (taus_, taus); |
428 |
|
|
429 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mets") != allNecessaryObjects.end()) |
429 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "mets") != objectsToGet.end()) |
430 |
|
event.getByLabel (mets_, mets); |
431 |
|
|
432 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "tracks") != allNecessaryObjects.end()) |
432 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "tracks") != objectsToGet.end()) |
433 |
|
event.getByLabel (tracks_, tracks); |
434 |
|
|
435 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "genjets") != allNecessaryObjects.end()) |
435 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "genjets") != objectsToGet.end()) |
436 |
|
event.getByLabel (genjets_, genjets); |
437 |
|
|
438 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "mcparticles") != allNecessaryObjects.end()) |
439 |
< |
event.getByLabel (mcparticles_, mcparticles); |
438 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "mcparticles") != objectsToGet.end()) |
439 |
> |
event.getByLabel (mcparticles_, mcparticles); |
440 |
> |
|
441 |
|
|
442 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "primaryvertexs") != allNecessaryObjects.end()) |
442 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "primaryvertexs") != objectsToGet.end()) |
443 |
|
event.getByLabel (primaryvertexs_, primaryvertexs); |
444 |
|
|
445 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "bxlumis") != allNecessaryObjects.end()) |
445 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "bxlumis") != objectsToGet.end()) |
446 |
|
event.getByLabel (bxlumis_, bxlumis); |
447 |
|
|
448 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "photons") != allNecessaryObjects.end()) |
448 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "photons") != objectsToGet.end()) |
449 |
|
event.getByLabel (photons_, photons); |
450 |
|
|
451 |
< |
if (std::find(allNecessaryObjects.begin(), allNecessaryObjects.end(), "superclusters") != allNecessaryObjects.end()) |
451 |
> |
if (std::find(objectsToGet.begin(), objectsToGet.end(), "superclusters") != objectsToGet.end()) |
452 |
|
event.getByLabel (superclusters_, superclusters); |
453 |
|
|
454 |
|
|
455 |
|
//get pile-up event weight |
456 |
|
double puScaleFactor = 1.00; |
457 |
< |
if(datasetType_ != "data") |
458 |
< |
puScaleFactor = puWeight_->at (events->at (0).numTruePV); |
457 |
> |
// if(datasetType_ != "data"){ |
458 |
> |
// puScaleFactor = puWeight_->at (events->at (0).numTruePV); |
459 |
> |
// cout << puScaleFactor << endl; |
460 |
> |
//} |
461 |
> |
|
462 |
> |
|
463 |
|
|
464 |
|
//loop over all channels |
465 |
|
|
477 |
|
} |
478 |
|
|
479 |
|
//loop over all cuts |
480 |
+ |
|
481 |
+ |
|
482 |
+ |
|
483 |
|
for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){ |
484 |
|
cut currentCut = currentChannel.cuts.at(currentCutIndex); |
485 |
|
|
486 |
< |
for(uint currentObjectIndex = 0; currentObjectIndex != allNecessaryObjects.size(); currentObjectIndex++){ |
487 |
< |
string currentObject = allNecessaryObjects.at(currentObjectIndex); |
486 |
> |
for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){ |
487 |
> |
string currentObject = objectsToCut.at(currentObjectIndex); |
488 |
|
|
489 |
+ |
//initialize maps to get ready to set cuts |
490 |
|
individualFlags[currentObject].push_back (vector<bool> ()); |
491 |
|
cumulativeFlags[currentObject].push_back (vector<bool> ()); |
492 |
|
|
493 |
+ |
} |
494 |
+ |
for(uint currentObjectIndex = 0; currentObjectIndex != objectsToCut.size(); currentObjectIndex++){ |
495 |
+ |
string currentObject = objectsToCut.at(currentObjectIndex); |
496 |
+ |
|
497 |
+ |
int flagsForPairCutsIndex = max(int(currentCutIndex-1),0); |
498 |
|
|
499 |
|
if(currentObject == "jets") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,jets.product(),"jets"); |
500 |
|
else if(currentObject == "muons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),"muons"); |
510 |
|
else if(currentObject == "photons") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,photons.product(),"photons"); |
511 |
|
else if(currentObject == "superclusters") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,superclusters.product(),"superclusters"); |
512 |
|
|
513 |
< |
else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(),"muon-muon pairs"); |
514 |
< |
else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(),"electron-electron pairs"); |
515 |
< |
else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(),"electron-muon pairs"); |
513 |
> |
|
514 |
> |
else if(currentObject == "muon-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,muons.product(),muons.product(), \ |
515 |
> |
cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \ |
516 |
> |
cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \ |
517 |
> |
"muon-muon pairs"); |
518 |
> |
else if(currentObject == "electron-electron pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),electrons.product(), \ |
519 |
> |
cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \ |
520 |
> |
cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \ |
521 |
> |
"electron-electron pairs"); |
522 |
> |
else if(currentObject == "electron-muon pairs") setObjectFlags(currentCut,currentCutIndex,individualFlags,cumulativeFlags,electrons.product(),muons.product(), \ |
523 |
> |
cumulativeFlags.at("electrons").at(flagsForPairCutsIndex), \ |
524 |
> |
cumulativeFlags.at("muons").at(flagsForPairCutsIndex), \ |
525 |
> |
"electron-muon pairs"); |
526 |
> |
|
527 |
|
|
528 |
|
} |
529 |
|
|
532 |
|
}//end loop over all cuts |
533 |
|
|
534 |
|
|
494 |
– |
|
535 |
|
//use cumulative flags to apply cuts at event level |
536 |
|
|
537 |
|
bool eventPassedAllCuts = true; |
540 |
|
eventPassedAllCuts = eventPassedAllCuts && triggerDecision; |
541 |
|
|
542 |
|
|
543 |
+ |
|
544 |
|
for(uint currentCutIndex = 0; currentCutIndex != currentChannel.cuts.size(); currentCutIndex++){ |
545 |
|
|
546 |
|
//loop over all objects and count how many passed the cumulative selection up to this point |
565 |
|
if(!eventPassedAllCuts)continue; |
566 |
|
|
567 |
|
|
568 |
< |
|
568 |
> |
|
569 |
|
|
570 |
|
//set position of primary vertex in event, in order to calculate quantities relative to it |
571 |
< |
vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back(); |
572 |
< |
for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){ |
573 |
< |
if(!vertexFlags.at(vertexIndex)) continue; |
574 |
< |
chosenPrimaryVertex = & primaryvertexs->at(vertexIndex); |
575 |
< |
break; |
571 |
> |
if(std::find(objectsToCut.begin(), objectsToCut.end(), "primaryvertexs") != objectsToCut.end()) { |
572 |
> |
vector<bool> vertexFlags = cumulativeFlags.at("primaryvertexs").back(); |
573 |
> |
for (uint vertexIndex = 0; vertexIndex != vertexFlags.size(); vertexIndex++){ |
574 |
> |
if(!vertexFlags.at(vertexIndex)) continue; |
575 |
> |
chosenPrimaryVertex = & primaryvertexs->at(vertexIndex); |
576 |
> |
break; |
577 |
> |
} |
578 |
> |
} |
579 |
> |
else { |
580 |
> |
chosenPrimaryVertex = & primaryvertexs->at(0); |
581 |
|
} |
582 |
|
|
583 |
+ |
|
584 |
+ |
|
585 |
|
//filling histograms |
586 |
|
for (uint histogramIndex = 0; histogramIndex != histograms.size(); histogramIndex++){ |
587 |
|
histogram currentHistogram = histograms.at(histogramIndex); |
589 |
|
if(currentHistogram.inputVariables.size() == 1){ |
590 |
|
TH1D* histo; |
591 |
|
histo = oneDHists_.at(currentChannelIndex).at(currentHistogram.name); |
592 |
< |
// std::cout << currentHistogram.inputCollection << std::endl; |
593 |
< |
// std::cout << "\t" << currentHistogram.inputVariables.at(0) << std::endl; |
592 |
> |
|
593 |
> |
|
594 |
> |
|
595 |
|
if(currentHistogram.inputCollection == "jets") fill1DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor); |
596 |
|
else if(currentHistogram.inputCollection == "muons") fill1DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor); |
597 |
< |
else if(currentHistogram.inputCollection == "muon-muon pairs") fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(),\ |
598 |
< |
cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(),\ |
599 |
< |
cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor); |
597 |
> |
else if(currentHistogram.inputCollection == "muon-muon pairs") |
598 |
> |
{ |
599 |
> |
fill1DHistogram(histo,currentHistogram,muons.product(),muons.product(), \ |
600 |
> |
cumulativeFlags.at("muons").back(),cumulativeFlags.at("muons").back(), \ |
601 |
> |
|
602 |
> |
cumulativeFlags.at("muon-muon pairs").back(),puScaleFactor); |
603 |
> |
|
604 |
> |
} |
605 |
|
else if(currentHistogram.inputCollection == "electrons") fill1DHistogram(histo,currentHistogram,electrons.product(),cumulativeFlags.at("electrons").back(),puScaleFactor); |
606 |
|
else if(currentHistogram.inputCollection == "electron-electron pairs") fill1DHistogram(histo,currentHistogram,electrons.product(),electrons.product(),\ |
607 |
|
cumulativeFlags.at("electrons").back(),cumulativeFlags.at("electrons").back(),\ |
623 |
|
else if(currentHistogram.inputVariables.size() == 2){ |
624 |
|
TH2D* histo; |
625 |
|
histo = twoDHists_.at(currentChannelIndex).at(currentHistogram.name); |
626 |
+ |
|
627 |
+ |
|
628 |
+ |
|
629 |
|
if(currentHistogram.inputCollection == "jets") fill2DHistogram(histo,currentHistogram,jets.product(),cumulativeFlags.at("jets").back(),puScaleFactor); |
630 |
|
else if(currentHistogram.inputCollection == "muons") fill2DHistogram(histo,currentHistogram,muons.product(),cumulativeFlags.at("muons").back(),puScaleFactor); |
631 |
|
else if(currentHistogram.inputCollection == "muon-muon pairs") fill2DHistogram(histo,currentHistogram,muons.product(),muons.product(), \ |
655 |
|
for (uint currentObjectIndex = 0; currentObjectIndex != objectsToPlot.size(); currentObjectIndex++){ |
656 |
|
string currentObject = objectsToPlot.at(currentObjectIndex); |
657 |
|
|
658 |
< |
if(currentObject == "muon-muon pairs") currentObject = "dimuonPairs"; |
659 |
< |
else if(currentObject == "electron-electron pairs") currentObject = "dielectronPairs"; |
660 |
< |
else if(currentObject == "electron-muon pairs") currentObject = "electronMuonPairs"; |
661 |
< |
string tempCurrentObject = currentObject; |
658 |
> |
string objectToPlot = ""; |
659 |
> |
|
660 |
> |
if(currentObject == "muon-muon pairs") objectToPlot = "dimuonPairs"; |
661 |
> |
else if(currentObject == "electron-electron pairs") objectToPlot = "dielectronPairs"; |
662 |
> |
else if(currentObject == "electron-muon pairs") objectToPlot = "electronMuonPairs"; |
663 |
> |
else objectToPlot = currentObject; |
664 |
> |
string tempCurrentObject = objectToPlot; |
665 |
|
tempCurrentObject.at(0) = toupper(tempCurrentObject.at(0)); |
666 |
|
string histoName = "num" + tempCurrentObject; |
667 |
|
|
668 |
|
|
669 |
< |
if(currentObject == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor); |
670 |
< |
else if(currentObject == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor); |
671 |
< |
else if(currentObject == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor); |
672 |
< |
else if(currentObject == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor); |
673 |
< |
else if(currentObject == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor); |
674 |
< |
else if(currentObject == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor); |
675 |
< |
else if(currentObject == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor); |
676 |
< |
else if(currentObject == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor); |
677 |
< |
else if(currentObject == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor); |
678 |
< |
else if(currentObject == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor); |
679 |
< |
else if(currentObject == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor); |
680 |
< |
else if(currentObject == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor); |
681 |
< |
else if(currentObject == "primaryvertexs"){ |
669 |
> |
|
670 |
> |
|
671 |
> |
//set position of primary vertex in event, in order to calculate quantities relative to it |
672 |
> |
if(std::find(objectsToCut.begin(), objectsToCut.end(), currentObject) != objectsToCut.end()) { |
673 |
> |
vector<bool> lastCutFlags = cumulativeFlags.at(currentObject).back(); |
674 |
> |
int numToPlot = 0; |
675 |
> |
for (uint currentFlag = 0; currentFlag != lastCutFlags.size(); currentFlag++){ |
676 |
> |
if(lastCutFlags.at(currentFlag)) numToPlot++; |
677 |
> |
} |
678 |
> |
if(objectToPlot == "primaryvertexs"){ |
679 |
> |
oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size()); |
680 |
> |
oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor); |
681 |
> |
} |
682 |
> |
else |
683 |
> |
oneDHists_.at(currentChannelIndex).at(histoName)->Fill(numToPlot,puScaleFactor); |
684 |
> |
} |
685 |
> |
else if(objectToPlot == "jets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(jets->size(),puScaleFactor); |
686 |
> |
else if(objectToPlot == "muons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size(),puScaleFactor); |
687 |
> |
else if(objectToPlot == "dimuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(muons->size()*(muons->size()-1)/2,puScaleFactor); |
688 |
> |
else if(objectToPlot == "electrons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size(),puScaleFactor); |
689 |
> |
else if(objectToPlot == "dielectronPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*(electrons->size()-1)/2,puScaleFactor); |
690 |
> |
else if(objectToPlot == "electronMuonPairs") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(electrons->size()*muons->size(),puScaleFactor); |
691 |
> |
else if(objectToPlot == "events") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(events->size(),puScaleFactor); |
692 |
> |
else if(objectToPlot == "taus") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(taus->size(),puScaleFactor); |
693 |
> |
else if(objectToPlot == "mets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mets->size(),puScaleFactor); |
694 |
> |
else if(objectToPlot == "tracks") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(tracks->size(),puScaleFactor); |
695 |
> |
else if(objectToPlot == "genjets") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(genjets->size(),puScaleFactor); |
696 |
> |
else if(objectToPlot == "mcparticles") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(mcparticles->size(),puScaleFactor); |
697 |
> |
else if(objectToPlot == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor); |
698 |
> |
else if(objectToPlot == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor); |
699 |
> |
else if(objectToPlot == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor); |
700 |
> |
else if(objectToPlot == "primaryvertexs"){ |
701 |
|
oneDHists_.at(currentChannelIndex).at(histoName+"BeforePileupCorrection")->Fill(primaryvertexs->size()); |
702 |
|
oneDHists_.at(currentChannelIndex).at(histoName+"AfterPileupCorrection")->Fill(primaryvertexs->size(),puScaleFactor); |
703 |
|
} |
625 |
– |
else if(currentObject == "bxlumis") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(bxlumis->size(),puScaleFactor); |
626 |
– |
else if(currentObject == "photons") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(photons->size(),puScaleFactor); |
627 |
– |
else if(currentObject == "superclusters") oneDHists_.at(currentChannelIndex).at(histoName)->Fill(superclusters->size(),puScaleFactor); |
704 |
|
|
705 |
|
} |
706 |
|
|
1052 |
|
else if(variable == "detIso") value = (object->trackIsoDR03) / object->pt; |
1053 |
|
else if(variable == "relPFdBetaIso") value = (object->pfIsoR04SumChargedHadronPt + max(0.0, object->pfIsoR04SumNeutralHadronEt + object->pfIsoR04SumPhotonEt - 0.5*object->pfIsoR04SumPUPt)) / object->pt; |
1054 |
|
else if(variable == "relPFrhoIso") value = ( object->chargedHadronIso + max(0.0, object->neutralHadronIso + object->photonIso - object->AEffDr03*object->rhoPrime) ) / object->pt; |
1055 |
+ |
|
1056 |
+ |
|
1057 |
+ |
|
1058 |
+ |
else if(variable == "correctedD0VertexInEBPlus"){ |
1059 |
+ |
if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex; |
1060 |
+ |
else value = -999; |
1061 |
+ |
} |
1062 |
+ |
else if(variable == "correctedD0VertexOutEBPlus"){ |
1063 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex; |
1064 |
+ |
else value = -999; |
1065 |
+ |
} |
1066 |
+ |
else if(variable == "correctedD0VertexEEPlus"){ |
1067 |
+ |
if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0Vertex; |
1068 |
+ |
else value = -999; |
1069 |
+ |
} |
1070 |
+ |
|
1071 |
+ |
else if(variable == "correctedD0BeamspotInEBPlus"){ |
1072 |
+ |
if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0; |
1073 |
+ |
else value = -999; |
1074 |
+ |
} |
1075 |
+ |
else if(variable == "correctedD0BeamspotOutEBPlus"){ |
1076 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0; |
1077 |
+ |
else value = -999; |
1078 |
+ |
} |
1079 |
+ |
else if(variable == "correctedD0BeamspotEEPlus"){ |
1080 |
+ |
if(fabs(object->eta) > 1.479 && object->eta > 0) value = object->correctedD0; |
1081 |
+ |
else value = -999; |
1082 |
+ |
} |
1083 |
+ |
|
1084 |
+ |
else if(variable == "correctedD0VertexInEBMinus"){ |
1085 |
+ |
if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex; |
1086 |
+ |
else value = -999; |
1087 |
+ |
} |
1088 |
+ |
else if(variable == "correctedD0VertexOutEBMinus"){ |
1089 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex; |
1090 |
+ |
else value = -999; |
1091 |
+ |
} |
1092 |
+ |
else if(variable == "correctedD0VertexEEMinus"){ |
1093 |
+ |
if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0Vertex; |
1094 |
+ |
else value = -999; |
1095 |
+ |
} |
1096 |
+ |
|
1097 |
+ |
else if(variable == "correctedD0BeamspotInEBMinus"){ |
1098 |
+ |
if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0; |
1099 |
+ |
else value = -999; |
1100 |
+ |
} |
1101 |
+ |
else if(variable == "correctedD0BeamspotOutEBMinus"){ |
1102 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0; |
1103 |
+ |
else value = -999; |
1104 |
+ |
} |
1105 |
+ |
else if(variable == "correctedD0BeamspotEEMinus"){ |
1106 |
+ |
if(fabs(object->eta) > 1.479 && object->eta < 0) value = object->correctedD0; |
1107 |
+ |
else value = -999; |
1108 |
+ |
} |
1109 |
+ |
|
1110 |
+ |
|
1111 |
+ |
else if(variable == "correctedD0VertexInEBPositiveCharge"){ |
1112 |
+ |
if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex; |
1113 |
+ |
else value = -999; |
1114 |
+ |
} |
1115 |
+ |
else if(variable == "correctedD0VertexOutEBPositiveCharge"){ |
1116 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex; |
1117 |
+ |
else value = -999; |
1118 |
+ |
} |
1119 |
+ |
else if(variable == "correctedD0VertexEEPositiveCharge"){ |
1120 |
+ |
if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0Vertex; |
1121 |
+ |
else value = -999; |
1122 |
+ |
} |
1123 |
+ |
|
1124 |
+ |
else if(variable == "correctedD0BeamspotInEBPositiveCharge"){ |
1125 |
+ |
if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0; |
1126 |
+ |
else value = -999; |
1127 |
+ |
} |
1128 |
+ |
else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){ |
1129 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0; |
1130 |
+ |
else value = -999; |
1131 |
+ |
} |
1132 |
+ |
else if(variable == "correctedD0BeamspotEEPositiveCharge"){ |
1133 |
+ |
if(fabs(object->eta) > 1.479 && object->charge > 0) value = object->correctedD0; |
1134 |
+ |
else value = -999; |
1135 |
+ |
} |
1136 |
+ |
|
1137 |
+ |
else if(variable == "correctedD0VertexInEBNegativeCharge"){ |
1138 |
+ |
if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex; |
1139 |
+ |
else value = -999; |
1140 |
+ |
} |
1141 |
+ |
else if(variable == "correctedD0VertexOutEBNegativeCharge"){ |
1142 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex; |
1143 |
+ |
else value = -999; |
1144 |
+ |
} |
1145 |
+ |
else if(variable == "correctedD0VertexEENegativeCharge"){ |
1146 |
+ |
if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0Vertex; |
1147 |
+ |
else value = -999; |
1148 |
+ |
} |
1149 |
+ |
|
1150 |
+ |
else if(variable == "correctedD0BeamspotInEBNegativeCharge"){ |
1151 |
+ |
if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0; |
1152 |
+ |
else value = -999; |
1153 |
+ |
} |
1154 |
+ |
else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){ |
1155 |
+ |
if(fabs(object->eta) < 1.479 && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0; |
1156 |
+ |
else value = -999; |
1157 |
+ |
} |
1158 |
+ |
else if(variable == "correctedD0BeamspotEENegativeCharge"){ |
1159 |
+ |
if(fabs(object->eta) > 1.479 && object->charge < 0) value = object->correctedD0; |
1160 |
+ |
else value = -999; |
1161 |
+ |
} |
1162 |
+ |
|
1163 |
+ |
|
1164 |
|
else if(variable == "tightID") { |
1165 |
|
value = object->isGlobalMuon > 0 \ |
1166 |
|
&& object->isPFMuon > 0 \ |
1181 |
|
&& object->numberOfLayersWithMeasurement > 5; |
1182 |
|
} |
1183 |
|
|
1184 |
+ |
|
1185 |
+ |
|
1186 |
|
else if(variable == "genMatchedId"){ |
1187 |
|
int index = getGenMatchedParticleIndex(object); |
1188 |
|
if(index == -1) value = 0; |
1380 |
|
else if(variable == "correctedD0VertexSig") value = object->correctedD0Vertex / hypot (object->tkD0err, hypot (chosenPrimaryVertex->xError, chosenPrimaryVertex->yError)); |
1381 |
|
else if(variable == "detIso") value = (object->trackIso) / object->pt; |
1382 |
|
else if(variable == "relPFrhoIso") value = ( object->chargedHadronIsoDR03 + max(0.0, object->neutralHadronIsoDR03 + object->photonIsoDR03 - object->AEffDr03*object->rhoPrime) ) / object->pt; |
1383 |
< |
else if(variable == "correctedD0VertexInEB"){ |
1384 |
< |
if(fabs(object->eta) < 0.8) value = object->correctedD0Vertex; |
1383 |
> |
|
1384 |
> |
|
1385 |
> |
else if(variable == "correctedD0VertexEEPositiveChargeLowPt"){ |
1386 |
> |
if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt > 45) value = object->correctedD0Vertex; |
1387 |
> |
else value = -999; |
1388 |
> |
} |
1389 |
> |
else if(variable == "correctedD0VertexEEPositiveChargeHighPt"){ |
1390 |
> |
if(fabs(object->eta) > 1.479 && object->charge > 0 && object->pt < 45) value = object->correctedD0Vertex; |
1391 |
> |
else value = -999; |
1392 |
> |
} |
1393 |
> |
|
1394 |
> |
|
1395 |
> |
else if(variable == "correctedD0VertexInEBPlus"){ |
1396 |
> |
if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0Vertex; |
1397 |
> |
else value = -999; |
1398 |
> |
} |
1399 |
> |
else if(variable == "correctedD0VertexOutEBPlus"){ |
1400 |
> |
if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0Vertex; |
1401 |
> |
else value = -999; |
1402 |
> |
} |
1403 |
> |
else if(variable == "correctedD0VertexEEPlus"){ |
1404 |
> |
if(object->isEE && object->eta > 0) value = object->correctedD0Vertex; |
1405 |
> |
else value = -999; |
1406 |
> |
} |
1407 |
> |
|
1408 |
> |
else if(variable == "correctedD0BeamspotInEBPlus"){ |
1409 |
> |
if(fabs(object->eta) < 0.8 && object->eta > 0) value = object->correctedD0; |
1410 |
|
else value = -999; |
1411 |
|
} |
1412 |
< |
else if(variable == "correctedD0VertexOutEB"){ |
1413 |
< |
if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0Vertex; |
1412 |
> |
else if(variable == "correctedD0BeamspotOutEBPlus"){ |
1413 |
> |
if(object->isEB && fabs(object->eta) > 0.8 && object->eta > 0) value = object->correctedD0; |
1414 |
|
else value = -999; |
1415 |
|
} |
1416 |
< |
else if(variable == "correctedD0VertexEE"){ |
1417 |
< |
if(object->isEE) value = object->correctedD0Vertex; |
1416 |
> |
else if(variable == "correctedD0BeamspotEEPlus"){ |
1417 |
> |
if(object->isEE && object->eta > 0) value = object->correctedD0; |
1418 |
|
else value = -999; |
1419 |
|
} |
1420 |
|
|
1421 |
< |
else if(variable == "correctedD0BeamspotInEB"){ |
1422 |
< |
if(fabs(object->eta) < 0.8) value = object->correctedD0; |
1421 |
> |
else if(variable == "correctedD0VertexInEBMinus"){ |
1422 |
> |
if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0Vertex; |
1423 |
|
else value = -999; |
1424 |
|
} |
1425 |
< |
else if(variable == "correctedD0BeamspotOutEB"){ |
1426 |
< |
if(object->isEB && fabs(object->eta) > 0.8) value = object->correctedD0; |
1425 |
> |
else if(variable == "correctedD0VertexOutEBMinus"){ |
1426 |
> |
if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0Vertex; |
1427 |
|
else value = -999; |
1428 |
|
} |
1429 |
< |
else if(variable == "correctedD0BeamspotEE"){ |
1430 |
< |
if(object->isEE) value = object->correctedD0; |
1429 |
> |
else if(variable == "correctedD0VertexEEMinus"){ |
1430 |
> |
if(object->isEE && object->eta < 0) value = object->correctedD0Vertex; |
1431 |
|
else value = -999; |
1432 |
|
} |
1433 |
|
|
1434 |
< |
else if(variable == "correctedD0OriginInEB"){ |
1435 |
< |
if(fabs(object->eta) < 0.8) value = object->tkD0; |
1434 |
> |
else if(variable == "correctedD0BeamspotInEBMinus"){ |
1435 |
> |
if(fabs(object->eta) < 0.8 && object->eta < 0) value = object->correctedD0; |
1436 |
|
else value = -999; |
1437 |
|
} |
1438 |
< |
else if(variable == "correctedD0OriginOutEB"){ |
1439 |
< |
if(object->isEB && fabs(object->eta) > 0.8) value = object->tkD0; |
1438 |
> |
else if(variable == "correctedD0BeamspotOutEBMinus"){ |
1439 |
> |
if(object->isEB && fabs(object->eta) > 0.8 && object->eta < 0) value = object->correctedD0; |
1440 |
|
else value = -999; |
1441 |
|
} |
1442 |
< |
else if(variable == "correctedD0OriginEE"){ |
1443 |
< |
if(object->isEE) value = object->tkD0; |
1442 |
> |
else if(variable == "correctedD0BeamspotEEMinus"){ |
1443 |
> |
if(object->isEE && object->eta < 0) value = object->correctedD0; |
1444 |
|
else value = -999; |
1445 |
|
} |
1446 |
|
|
1447 |
+ |
|
1448 |
+ |
else if(variable == "correctedD0VertexInEBPositiveCharge"){ |
1449 |
+ |
if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0Vertex; |
1450 |
+ |
else value = -999; |
1451 |
+ |
} |
1452 |
+ |
else if(variable == "correctedD0VertexOutEBPositiveCharge"){ |
1453 |
+ |
if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0Vertex; |
1454 |
+ |
else value = -999; |
1455 |
+ |
} |
1456 |
+ |
else if(variable == "correctedD0VertexEEPositiveCharge"){ |
1457 |
+ |
if(object->isEE && object->charge > 0) value = object->correctedD0Vertex; |
1458 |
+ |
else value = -999; |
1459 |
+ |
} |
1460 |
+ |
|
1461 |
+ |
else if(variable == "correctedD0BeamspotInEBPositiveCharge"){ |
1462 |
+ |
if(fabs(object->eta) < 0.8 && object->charge > 0) value = object->correctedD0; |
1463 |
+ |
else value = -999; |
1464 |
+ |
} |
1465 |
+ |
else if(variable == "correctedD0BeamspotOutEBPositiveCharge"){ |
1466 |
+ |
if(object->isEB && fabs(object->eta) > 0.8 && object->charge > 0) value = object->correctedD0; |
1467 |
+ |
else value = -999; |
1468 |
+ |
} |
1469 |
+ |
else if(variable == "correctedD0BeamspotEEPositiveCharge"){ |
1470 |
+ |
if(object->isEE && object->charge > 0) value = object->correctedD0; |
1471 |
+ |
else value = -999; |
1472 |
+ |
} |
1473 |
+ |
|
1474 |
+ |
else if(variable == "correctedD0VertexInEBNegativeCharge"){ |
1475 |
+ |
if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0Vertex; |
1476 |
+ |
else value = -999; |
1477 |
+ |
} |
1478 |
+ |
else if(variable == "correctedD0VertexOutEBNegativeCharge"){ |
1479 |
+ |
if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0Vertex; |
1480 |
+ |
else value = -999; |
1481 |
+ |
} |
1482 |
+ |
else if(variable == "correctedD0VertexEENegativeCharge"){ |
1483 |
+ |
if(object->isEE && object->charge < 0) value = object->correctedD0Vertex; |
1484 |
+ |
else value = -999; |
1485 |
+ |
} |
1486 |
+ |
|
1487 |
+ |
else if(variable == "correctedD0BeamspotInEBNegativeCharge"){ |
1488 |
+ |
if(fabs(object->eta) < 0.8 && object->charge < 0) value = object->correctedD0; |
1489 |
+ |
else value = -999; |
1490 |
+ |
} |
1491 |
+ |
else if(variable == "correctedD0BeamspotOutEBNegativeCharge"){ |
1492 |
+ |
if(object->isEB && fabs(object->eta) > 0.8 && object->charge < 0) value = object->correctedD0; |
1493 |
+ |
else value = -999; |
1494 |
+ |
} |
1495 |
+ |
else if(variable == "correctedD0BeamspotEENegativeCharge"){ |
1496 |
+ |
if(object->isEE && object->charge < 0) value = object->correctedD0; |
1497 |
+ |
else value = -999; |
1498 |
+ |
} |
1499 |
+ |
|
1500 |
+ |
|
1501 |
|
else if(variable == "tightIDdisplaced"){ |
1502 |
|
if (object->isEB) |
1503 |
|
{ |
1776 |
|
OSUAnalysis::valueLookup (const BNtrack* object, string variable, string function){ |
1777 |
|
|
1778 |
|
double value = 0.0; |
1779 |
< |
|
1779 |
> |
double pMag = sqrt(object->pt * object->pt + |
1780 |
> |
object->pz * object->pz); |
1781 |
> |
|
1782 |
|
if(variable == "pt") value = object->pt; |
1783 |
|
else if(variable == "px") value = object->px; |
1784 |
|
else if(variable == "py") value = object->py; |
1798 |
|
else if(variable == "isHighPurity") value = object->isHighPurity; |
1799 |
|
|
1800 |
|
|
1801 |
+ |
//additional BNs info for disappTrks |
1802 |
+ |
else if(variable == "isGoodPtResolution") value = object->isGoodPtResolution; |
1803 |
+ |
else if(variable == "caloEMDeltaRp3") value = object->caloEMDeltaRp3; |
1804 |
+ |
else if(variable == "caloHadDeltaRp3") value = object->caloHadDeltaRp3; |
1805 |
+ |
else if(variable == "caloEMDeltaRp4") value = object->caloEMDeltaRp4; |
1806 |
+ |
else if(variable == "caloHadDeltaRp4") value = object->caloHadDeltaRp4; |
1807 |
+ |
else if(variable == "caloEMDeltaRp5") value = object->caloEMDeltaRp5; |
1808 |
+ |
else if(variable == "caloHadDeltaRp5") value = object->caloHadDeltaRp5; |
1809 |
+ |
else if(variable == "nHitsMissingOuter") value = object->nHitsMissingOuter; |
1810 |
+ |
else if(variable == "nHitsMissingInner") value = object->nHitsMissingInner; |
1811 |
+ |
else if(variable == "nHitsMissingMiddle") value = object->nHitsMissingMiddle; |
1812 |
+ |
//user defined variables |
1813 |
+ |
else if(variable == "d0wrtBS") value = (object->vx-events->at(0).BSx)*object->py/object->pt - (object->vy-events->at(0).BSy)*object->px/object->pt; |
1814 |
+ |
else if(variable == "dZwrtBS") value = object->dZ - events->at(0).BSz; |
1815 |
+ |
else if(variable == "caloTotDeltaRp5") value =(object->caloHadDeltaRp5 + object->caloEMDeltaRp5); |
1816 |
+ |
else if(variable == "caloTotDeltaRp5ByP") value =( (object->caloHadDeltaRp5 + object->caloEMDeltaRp5)/pMag); |
1817 |
+ |
else if(variable == "isIso") value = getTrkIsIso(object, tracks.product()); |
1818 |
+ |
else if(variable == "isMatchedDeadEcal") value = getTrkIsMatchedDeadEcal(object); |
1819 |
+ |
else if(variable == "ptErrorByPt") value = (object->ptError/object->pt); |
1820 |
+ |
else if(variable == "ptError") value = object->ptError; |
1821 |
+ |
else if(variable == "ptRes") value = getTrkPtRes(object); |
1822 |
+ |
|
1823 |
+ |
|
1824 |
|
else if(variable == "genMatchedId"){ |
1825 |
|
int index = getGenMatchedParticleIndex(object); |
1826 |
|
if(index == -1) value = 0; |
2189 |
|
|
2190 |
|
double value = 0.0; |
2191 |
|
|
2192 |
< |
if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi); |
2192 |
> |
if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi)); |
2193 |
> |
else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi); |
2194 |
|
else if(variable == "invMass"){ |
2195 |
|
TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy); |
2196 |
|
TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy); |
2197 |
|
|
2198 |
|
value = (fourVector1 + fourVector2).M();} |
2199 |
+ |
|
2200 |
|
else if(variable == "threeDAngle") |
2201 |
|
{ |
2202 |
|
TVector3 threeVector1(object1->px, object1->py, object1->pz); |
2203 |
|
TVector3 threeVector2(object2->px, object2->py, object2->pz); |
2204 |
|
value = (threeVector1.Angle(threeVector2)); |
2205 |
|
} |
2206 |
+ |
else if(variable == "alpha") |
2207 |
+ |
{ |
2208 |
+ |
static const double pi = 3.1415926535897932384626433832795028841971693993751058; |
2209 |
+ |
TVector3 threeVector1(object1->px, object1->py, object1->pz); |
2210 |
+ |
TVector3 threeVector2(object2->px, object2->py, object2->pz); |
2211 |
+ |
value = (pi-threeVector1.Angle( threeVector2)); |
2212 |
+ |
} |
2213 |
|
|
2214 |
+ |
else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi); |
2215 |
|
|
2216 |
|
|
2217 |
|
else if(variable == "deltaCorrectedD0Vertex") value = object1->correctedD0Vertex - object2->correctedD0Vertex; |
2228 |
|
else if(variable == "muon2CorrectedD0Vertex"){ |
2229 |
|
value = object2->correctedD0Vertex; |
2230 |
|
} |
2231 |
+ |
else if(variable == "muon1timeAtIpInOut"){ |
2232 |
+ |
value = object1->timeAtIpInOut; |
2233 |
+ |
} |
2234 |
+ |
else if(variable == "muon2timeAtIpInOut"){ |
2235 |
+ |
value = object2->timeAtIpInOut; |
2236 |
+ |
} |
2237 |
|
else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;} |
2238 |
|
|
2239 |
|
|
2253 |
|
|
2254 |
|
double value = 0.0; |
2255 |
|
|
2256 |
< |
if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi); |
2256 |
> |
if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi)); |
2257 |
> |
else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi); |
2258 |
|
else if(variable == "invMass"){ |
2259 |
|
TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy); |
2260 |
|
TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy); |
2302 |
|
|
2303 |
|
double value = 0.0; |
2304 |
|
|
2305 |
< |
if(variable == "deltaPhi") value = deltaPhi(object1->phi,object2->phi); |
2305 |
> |
if(variable == "deltaPhi") value = fabs(deltaPhi(object1->phi,object2->phi)); |
2306 |
> |
else if(variable == "deltaR") value = deltaR(object1->eta,object1->phi,object2->eta,object2->phi); |
2307 |
|
else if(variable == "invMass"){ |
2308 |
|
TLorentzVector fourVector1(object1->px, object1->py, object1->pz, object1->energy); |
2309 |
|
TLorentzVector fourVector2(object2->px, object2->py, object2->pz, object2->energy); |
2347 |
|
} |
2348 |
|
|
2349 |
|
|
2350 |
+ |
// Calculate the number of tracks in cone of DeltaR<0.5 around track1. |
2351 |
+ |
// Return true iff no other tracks are found in this cone. |
2352 |
+ |
int |
2353 |
+ |
OSUAnalysis::getTrkIsIso (const BNtrack* track1, const BNtrackCollection* trackColl){ |
2354 |
+ |
for(BNtrackCollection::const_iterator track2 = trackColl->begin(); track2 !=trackColl->end(); track2++){ |
2355 |
+ |
if(track1->eta == track2->eta && track1->phi == track2->phi) continue; // Do not compare the track to itself. |
2356 |
+ |
double deltaRtrk = deltaR(track1->eta, track1->phi, track2->eta, track2->phi); |
2357 |
+ |
if(deltaRtrk < 0.5) return 0; |
2358 |
+ |
} |
2359 |
+ |
return 1; |
2360 |
+ |
|
2361 |
+ |
} |
2362 |
+ |
|
2363 |
+ |
|
2364 |
+ |
double |
2365 |
+ |
OSUAnalysis::getTrkPtRes (const BNtrack* track1){ |
2366 |
+ |
|
2367 |
+ |
double ptTrue = getTrkPtTrue(track1, mcparticles.product()); |
2368 |
+ |
double PtRes = (track1->pt - ptTrue) / ptTrue; |
2369 |
+ |
|
2370 |
+ |
return PtRes; |
2371 |
+ |
|
2372 |
+ |
} |
2373 |
+ |
|
2374 |
+ |
|
2375 |
+ |
double |
2376 |
+ |
OSUAnalysis::getTrkPtTrue (const BNtrack* track1, const BNmcparticleCollection* genPartColl){ |
2377 |
+ |
double value = -99; |
2378 |
+ |
double genDeltaRLowest = 999; |
2379 |
+ |
|
2380 |
+ |
for (BNmcparticleCollection::const_iterator genPart = genPartColl->begin(); genPart !=genPartColl->end(); genPart++){ |
2381 |
+ |
double genDeltaRtemp = deltaR(genPart->eta, genPart->phi,track1->eta, track1->phi); |
2382 |
+ |
if (genDeltaRtemp < genDeltaRLowest) { |
2383 |
+ |
genDeltaRLowest = genDeltaRtemp; |
2384 |
+ |
if (genDeltaRLowest < 0.05) { // Only consider it truth-matched if DeltaR<0.15. |
2385 |
+ |
double ptTrue = genPart->pt; |
2386 |
+ |
value = ptTrue; |
2387 |
+ |
} |
2388 |
+ |
} |
2389 |
+ |
} |
2390 |
+ |
|
2391 |
+ |
return value; |
2392 |
+ |
|
2393 |
+ |
} |
2394 |
+ |
|
2395 |
+ |
//creates a map of the dead Ecal channels in the barrel and endcap |
2396 |
+ |
//to see how the map of dead Ecal channels is created look at function getChannelStatusMaps() here: |
2397 |
+ |
//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/jbrinson/DisappTrk/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc?revision=1.88&view=markup |
2398 |
+ |
void |
2399 |
+ |
OSUAnalysis::WriteDeadEcal (){ |
2400 |
+ |
double etaEcal, phiEcal; |
2401 |
+ |
ifstream DeadEcalFile(deadEcalFile_); |
2402 |
+ |
if(!DeadEcalFile) { |
2403 |
+ |
cout << "Error: DeadEcalFile has not been found." << endl; |
2404 |
+ |
return; |
2405 |
+ |
} |
2406 |
+ |
if(DeadEcalVec.size()!= 0){ |
2407 |
+ |
cout << "Error: DeadEcalVec has a nonzero size" << endl; |
2408 |
+ |
return; |
2409 |
+ |
} |
2410 |
+ |
while(!DeadEcalFile.eof()) |
2411 |
+ |
{ |
2412 |
+ |
DeadEcalFile >> etaEcal >> phiEcal; |
2413 |
+ |
DeadEcal newChan; |
2414 |
+ |
newChan.etaEcal = etaEcal; |
2415 |
+ |
newChan.phiEcal = phiEcal; |
2416 |
+ |
DeadEcalVec.push_back(newChan); |
2417 |
+ |
} |
2418 |
+ |
if(DeadEcalVec.size() == 0) cout << "Warning: No dead Ecal channels have been found." << endl; |
2419 |
+ |
} |
2420 |
+ |
|
2421 |
+ |
//if a track is found within dR<0.05 of a dead Ecal channel value = 1, otherwise value = 0 |
2422 |
+ |
int |
2423 |
+ |
OSUAnalysis::getTrkIsMatchedDeadEcal (const BNtrack* track1){ |
2424 |
+ |
double deltaRLowest = 999; |
2425 |
+ |
int value = 0; |
2426 |
+ |
if (DeadEcalVec.size() == 0) WriteDeadEcal(); |
2427 |
+ |
for(std::vector<DeadEcal>::const_iterator ecal = DeadEcalVec.begin(); ecal != DeadEcalVec.end(); ++ecal){ |
2428 |
+ |
double eta = ecal->etaEcal; |
2429 |
+ |
double phi = ecal->phiEcal; |
2430 |
+ |
double deltaRtemp = deltaR(eta, track1->eta, phi, track1->phi); |
2431 |
+ |
if(deltaRtemp < deltaRLowest) deltaRLowest = deltaRtemp; |
2432 |
+ |
} |
2433 |
+ |
if (deltaRLowest<0.05) {value = 1;} |
2434 |
+ |
else {value = 0;} |
2435 |
+ |
return value; |
2436 |
+ |
} |
2437 |
+ |
|
2438 |
+ |
|
2439 |
+ |
|
2440 |
+ |
|
2441 |
|
double |
2442 |
|
OSUAnalysis::applyFunction(string function, double value){ |
2443 |
|
|
2444 |
|
if(function == "abs") value = fabs(value); |
2445 |
|
else if(function == "fabs") value = fabs(value); |
2446 |
< |
|
2446 |
> |
else if(function == "log10") value = log10(value); |
2447 |
> |
else if(function == "log") value = log10(value); |
2448 |
|
|
2449 |
|
else if(function == "") value = value; |
2450 |
|
else{std::cout << "WARNING: invalid function '" << function << "'\n";} |
2457 |
|
template <class InputCollection> |
2458 |
|
void OSUAnalysis::setObjectFlags(cut ¤tCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, InputCollection inputCollection, string inputType){ |
2459 |
|
|
2460 |
+ |
|
2461 |
|
for (uint object = 0; object != inputCollection->size(); object++){ |
2462 |
|
|
2463 |
+ |
|
2464 |
|
bool decision = true;//object passes if this cut doesn't cut on that type of object |
2465 |
|
|
2466 |
+ |
|
2467 |
|
if(currentCut.inputCollection == inputType){ |
2468 |
|
|
2469 |
|
vector<bool> subcutDecisions; |
2485 |
|
} |
2486 |
|
individualFlags.at(inputType).at(currentCutIndex).push_back(decision); |
2487 |
|
|
2488 |
+ |
|
2489 |
|
//set flags for objects that pass each cut AND all the previous cuts |
2490 |
|
bool previousCumulativeFlag = true; |
2491 |
|
for(uint previousCutIndex = 0; previousCutIndex != currentCutIndex; previousCutIndex++){ |
2494 |
|
} |
2495 |
|
cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision); |
2496 |
|
|
2497 |
+ |
|
2498 |
|
} |
2499 |
|
|
2500 |
|
} |
2502 |
|
|
2503 |
|
template <class InputCollection1, class InputCollection2> |
2504 |
|
void OSUAnalysis::setObjectFlags(cut ¤tCut, uint currentCutIndex, flagMap &individualFlags, flagMap &cumulativeFlags, \ |
2505 |
< |
InputCollection1 inputCollection1, InputCollection2 inputCollection2, string inputType){ |
2505 |
> |
InputCollection1 inputCollection1, InputCollection2 inputCollection2, vector<bool> flags1, vector<bool> flags2, string inputType){ |
2506 |
|
|
2507 |
|
|
2508 |
|
bool sameObjects = false; |
2509 |
|
if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true; |
2510 |
|
|
2511 |
+ |
|
2512 |
|
int counter = 0; |
2513 |
|
for (uint object1 = 0; object1 != inputCollection1->size(); object1++){ |
2514 |
|
for (uint object2 = 0; object2 != inputCollection2->size(); object2++){ |
2546 |
|
if(previousCumulativeFlag && individualFlags.at(inputType).at(previousCutIndex).at(counter)) previousCumulativeFlag = true; |
2547 |
|
else{ previousCumulativeFlag = false; break;} |
2548 |
|
} |
2549 |
< |
cumulativeFlags.at(inputType).at(currentCutIndex).push_back(previousCumulativeFlag && decision); |
2549 |
> |
//apply flags for the components of the composite object as well |
2550 |
> |
bool currentCumulativeFlag = true; |
2551 |
> |
if(flags1.size() == 0 && flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision; |
2552 |
> |
else if(flags1.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags2.at(object2); |
2553 |
> |
else if(flags2.size() == 0) currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1); |
2554 |
> |
else currentCumulativeFlag = previousCumulativeFlag && decision && flags1.at(object1) && flags2.at(object2); |
2555 |
> |
cumulativeFlags.at(inputType).at(currentCutIndex).push_back(currentCumulativeFlag); |
2556 |
|
|
2557 |
|
counter++; |
2558 |
|
} |
2595 |
|
bool sameObjects = false; |
2596 |
|
if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true; |
2597 |
|
|
2598 |
< |
int counter = 0; |
2598 |
> |
int pairCounter = 0; |
2599 |
|
for (uint object1 = 0; object1 != inputCollection1->size(); object1++){ |
2600 |
|
for (uint object2 = 0; object2 != inputCollection2->size(); object2++){ |
2601 |
|
|
2604 |
|
//only take objects which have passed all cuts and pairs which have passed all cuts |
2605 |
|
if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue; |
2606 |
|
if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue; |
2607 |
< |
if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue; |
2607 |
> |
if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue; |
2608 |
|
|
2609 |
|
string currentString = parameters.inputVariables.at(0); |
2610 |
|
string inputVariable = ""; |
2621 |
|
double value = valueLookup(&inputCollection1->at(object1), &inputCollection2->at(object2), inputVariable, function); |
2622 |
|
histo->Fill(value,puScaleFactor); |
2623 |
|
|
2624 |
< |
counter++; |
2624 |
> |
pairCounter++; |
2625 |
|
} |
2626 |
|
} |
2627 |
|
|
2674 |
|
bool sameObjects = false; |
2675 |
|
if(typeid(InputCollection1).name() == typeid(InputCollection2).name()) sameObjects = true; |
2676 |
|
|
2677 |
< |
int counter = 0; |
2677 |
> |
int pairCounter = 0; |
2678 |
|
for (uint object1 = 0; object1 != inputCollection1->size(); object1++){ |
2679 |
|
for (uint object2 = 0; object2 != inputCollection2->size(); object2++){ |
2680 |
|
|
2683 |
|
//only take objects which have passed all cuts and pairs which have passed all cuts |
2684 |
|
if(!plotAllObjectsInPassingEvents_ && !flags1.at(object1)) continue; |
2685 |
|
if(!plotAllObjectsInPassingEvents_ && !flags2.at(object2)) continue; |
2686 |
< |
if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(counter)) continue; |
2686 |
> |
if(!plotAllObjectsInPassingEvents_ && !pairFlags.at(pairCounter)) continue; |
2687 |
|
|
2688 |
|
string currentString = parameters.inputVariables.at(0); |
2689 |
|
string inputVariable = ""; |
2714 |
|
|
2715 |
|
histo->Fill(valueX,valueY,puScaleFactor); |
2716 |
|
|
2717 |
< |
counter++; |
2717 |
> |
pairCounter++; |
2718 |
> |
|
2719 |
|
} |
2720 |
|
} |
2721 |
|
|
2724 |
|
|
2725 |
|
template <class InputObject> |
2726 |
|
int OSUAnalysis::getGenMatchedParticleIndex(InputObject object){ |
2727 |
< |
// std::cout << "starting gen match\n"; |
2727 |
> |
|
2728 |
|
int bestMatchIndex = -1; |
2729 |
|
double bestMatchDeltaR = 999; |
2316 |
– |
//std::cout << "\n\n\n"; |
2317 |
– |
for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){ |
2318 |
– |
//std::cout << "pdgId = " << mcparticle->id << " status = " << mcparticle->status << " motherId = " << mcparticle->motherId << " grandMotherId = " << mcparticle->grandMotherId << std::endl; |
2730 |
|
|
2731 |
+ |
for(BNmcparticleCollection::const_iterator mcparticle = mcparticles->begin (); mcparticle != mcparticles->end (); mcparticle++){ |
2732 |
|
|
2733 |
|
double currentDeltaR = deltaR(object->eta,object->phi,mcparticle->eta,mcparticle->phi); |
2734 |
|
if(currentDeltaR > 0.05) continue; |
2739 |
|
} |
2740 |
|
|
2741 |
|
} |
2742 |
< |
// if(bestMatchIndex != -1) std::cout << "best match: pdgId = " << mcparticles->at(bestMatchIndex).id << " status = " << mcparticles->at(bestMatchIndex).status << " motherId = " << mcparticles->at(bestMatchIndex).motherId << " grandMotherId = " << mcparticles->at(bestMatchIndex).grandMotherId << std::endl; |
2331 |
< |
// std::cout << "finished gen match\n"; |
2742 |
> |
|
2743 |
|
return bestMatchIndex; |
2744 |
|
|
2745 |
|
} |
2763 |
|
} |
2764 |
|
|
2765 |
|
} |
2766 |
< |
// if(bestMatchIndex != -1) std::cout << "best match: pdgId = " << mcparticles->at(bestMatchIndex).id << " status = " << mcparticles->at(bestMatchIndex).status << " motherId = " << mcparticles->at(bestMatchIndex).motherId << " grandMotherId = " << mcparticles->at(bestMatchIndex).grandMotherId << std::endl; |
2766 |
> |
|
2767 |
|
return bestMatchIndex; |
2768 |
|
} |
2769 |
|
|