17 #include "TClonesArray.h"
22 #include <boost/regex.hpp>
29 double NPV =
static_cast<double>(intNPV);
30 double etaArray[101] = {-5, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4, -3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5};
32 for(
int i=0;
i<100; ++
i){
33 if(eta > etaArray[
i] && eta < etaArray[
i+1]){
41 double p0[100] = {0.08187, 0.096718, 0.11565, 0.12115, 0.12511, 0.12554, 0.13858, 0.14282, 0.14302, 0.15054, 0.14136, 0.14992, 0.13812, 0.13771, 0.13165, 0.12609, 0.12446, 0.11311, 0.13771, 0.16401, 0.20454, 0.27899, 0.34242, 0.43096, 0.50742, 0.59683, 0.66877, 0.68664, 0.69541, 0.66873, 0.64175, 0.61097, 0.58524, 0.5712, 0.55752, 0.54869, 0.31073, 0.22667, 0.55614, 0.55962, 0.54348, 0.53206, 0.51594, 0.49928, 0.49139, 0.48766, 0.49157, 0.49587, 0.50109, 0.5058, 0.51279, 0.51515, 0.51849, 0.52607, 0.52362, 0.52169, 0.53579, 0.54821, 0.56262, 0.58355, 0.58809, 0.57525, 0.52539, 0.53505, 0.52307, 0.52616, 0.52678, 0.53536, 0.55141, 0.58107, 0.60556, 0.62601, 0.60897, 0.59018, 0.49593, 0.40462, 0.32052, 0.24436, 0.18867, 0.12591, 0.095421, 0.090578, 0.078767, 0.11797, 0.14057, 0.14614, 0.15232, 0.14742, 0.15647, 0.14947, 0.15805, 0.14467, 0.14526, 0.14081, 0.1262, 0.12429, 0.11951, 0.11146, 0.095677, 0.083126};
42 double p1[100] = {0.26831, 0.30901, 0.37017, 0.38747, 0.41547, 0.45237, 0.49963, 0.54074, 0.54949, 0.5937, 0.56904, 0.60766, 0.58042, 0.59823, 0.58535, 0.54594, 0.58403, 0.601, 0.65401, 0.65049, 0.65264, 0.6387, 0.60646, 0.59669, 0.55561, 0.5053, 0.42889, 0.37264, 0.36456, 0.36088, 0.36728, 0.37439, 0.38779, 0.40133, 0.40989, 0.41722, 0.47539, 0.49848, 0.42642, 0.42431, 0.42113, 0.41285, 0.41003, 0.41116, 0.41231, 0.41634, 0.41795, 0.41806, 0.41786, 0.41765, 0.41779, 0.41961, 0.42144, 0.42192, 0.4209, 0.41885, 0.4163, 0.4153, 0.41864, 0.4257, 0.43018, 0.43218, 0.43798, 0.42723, 0.42185, 0.41349, 0.40553, 0.39132, 0.3779, 0.37055, 0.36522, 0.37057, 0.38058, 0.43259, 0.51052, 0.55918, 0.60178, 0.60995, 0.64087, 0.65554, 0.65308, 0.65654, 0.60466, 0.58678, 0.54392, 0.58277, 0.59651, 0.57916, 0.60744, 0.56882, 0.59323, 0.5499, 0.54003, 0.49938, 0.4511, 0.41499, 0.38676, 0.36955, 0.30803, 0.26659};
43 double p2[100] = {0.00080918, 0.00083447, 0.0011378, 0.0011221, 0.0013613, 0.0016362, 0.0015854, 0.0019131, 0.0017474, 0.0020078, 0.001856, 0.0020331, 0.0020823, 0.001898, 0.0020096, 0.0016464, 0.0032413, 0.0045615, 0.0054495, 0.0057584, 0.0058982, 0.0058956, 0.0055109, 0.0051433, 0.0042098, 0.0032096, 0.00044089, -0.0003884, -0.0007059, -0.00092769, -0.001116, -0.0010437, -0.00080318, -0.00044142, 6.7232e-05, 0.00055265, -0.0014486, -0.0020432, 0.0015121, 0.0016343, 0.0015638, 0.0015707, 0.0014403, 0.0012886, 0.0011684, 0.00099089, 0.00091497, 0.00087915, 0.00084703, 0.00084542, 0.00087419, 0.00088013, 0.00090493, 0.00095853, 0.0010389, 0.0011191, 0.0012643, 0.0013833, 0.001474, 0.0015401, 0.0015582, 0.0014265, 0.00087453, 0.00086639, 0.00042986, -5.0257e-06, -0.00053124, -0.00086417, -0.0011228, -0.0011749, -0.0010068, -0.00083012, -0.00062906, 0.00021515, 0.0028714, 0.0038835, 0.0047212, 0.0051427, 0.0055762, 0.0055872, 0.0054989, 0.0053033, 0.0044519, 0.0032223, 0.0017641, 0.0021165, 0.0019909, 0.0021061, 0.0020322, 0.0018357, 0.0019829, 0.001683, 0.0018553, 0.0015304, 0.0015822, 0.0013119, 0.0010745, 0.0010808, 0.00080678, 0.00079756};
44 pt_density = p0[ind] + p1[ind]*(NPV - 1) + p2[ind]*(NPV - 1)*(NPV - 1);
47 double p0[100] = {0.12523, 0.14896, 0.17696, 0.19376, 0.20038, 0.21353, 0.25069, 0.27089, 0.29124, 0.31947, 0.31781, 0.35453, 0.35424, 0.38159, 0.39453, 0.4003, 0.34798, 0.26303, 0.24824, 0.22857, 0.22609, 0.26793, 0.30096, 0.37637, 0.44461, 0.55692, 0.70328, 0.72458, 0.75065, 0.73569, 0.72485, 0.69933, 0.69804, 0.70775, 0.70965, 0.71439, 0.72189, 0.73691, 0.74847, 0.74968, 0.73467, 0.70115, 0.6732, 0.65971, 0.65724, 0.67751, 0.69569, 0.70905, 0.71815, 0.72119, 0.72128, 0.71645, 0.70588, 0.68958, 0.66978, 0.65959, 0.66889, 0.68713, 0.71063, 0.74283, 0.75153, 0.74733, 0.73335, 0.71346, 0.70168, 0.69445, 0.68841, 0.67761, 0.67654, 0.6957, 0.70276, 0.71057, 0.68176, 0.64651, 0.49156, 0.38366, 0.31375, 0.24127, 0.21395, 0.17783, 0.19026, 0.21486, 0.24689, 0.3434, 0.40184, 0.39876, 0.3873, 0.36462, 0.36337, 0.32777, 0.328, 0.29868, 0.28087, 0.25713, 0.22466, 0.20784, 0.19798, 0.18054, 0.15022, 0.12811};
48 double p1[100] = {0.26829, 0.30825, 0.37034, 0.38736, 0.41645, 0.45985, 0.51433, 0.56215, 0.5805, 0.63926, 0.62007, 0.67895, 0.66015, 0.68817, 0.67975, 0.64161, 0.70887, 0.74454, 0.80197, 0.78873, 0.77892, 0.74943, 0.70034, 0.6735, 0.60774, 0.53312, 0.42132, 0.36279, 0.3547, 0.35014, 0.35655, 0.3646, 0.37809, 0.38922, 0.39599, 0.40116, 0.40468, 0.40645, 0.40569, 0.4036, 0.39874, 0.39326, 0.39352, 0.39761, 0.40232, 0.40729, 0.41091, 0.41247, 0.413, 0.41283, 0.41289, 0.4134, 0.41322, 0.41185, 0.40769, 0.40193, 0.39707, 0.39254, 0.39274, 0.3989, 0.40474, 0.40758, 0.40788, 0.40667, 0.40433, 0.40013, 0.39371, 0.38154, 0.36723, 0.3583, 0.35148, 0.35556, 0.36172, 0.41073, 0.50629, 0.57068, 0.62972, 0.65188, 0.69954, 0.72967, 0.74333, 0.76148, 0.71418, 0.69062, 0.63065, 0.67117, 0.68278, 0.66028, 0.68147, 0.62494, 0.64452, 0.58685, 0.57076, 0.52387, 0.47132, 0.42637, 0.39554, 0.37989, 0.31825, 0.27969};
49 double p2[100] = {-0.0014595, -0.0014618, -0.0011988, -0.00095404, -5.3893e-05, 0.00018901, 0.00012553, 0.0004172, 0.00020229, 0.00051942, 0.00052088, 0.00076727, 0.0010407, 0.0010184, 0.0013442, 0.0011271, 0.0032841, 0.0045259, 0.0051803, 0.0054477, 0.0055691, 0.0056668, 0.0053084, 0.0050978, 0.0042061, 0.003321, 0.00045155, 0.00021376, 0.0001178, -2.6836e-05, -0.00017689, -0.00014723, 0.00016887, 0.00067322, 0.0012952, 0.0019229, 0.0024702, 0.0028854, 0.0031576, 0.003284, 0.0032643, 0.0031061, 0.0028377, 0.0025386, 0.0022583, 0.0020448, 0.001888, 0.0017968, 0.0017286, 0.0016989, 0.0017014, 0.0017302, 0.0017958, 0.0018891, 0.0020609, 0.0022876, 0.0025391, 0.0028109, 0.0030294, 0.0031867, 0.0032068, 0.0030755, 0.0028181, 0.0023893, 0.0018359, 0.0012192, 0.00061654, 0.00016088, -0.00015204, -0.00019503, -3.7236e-05, 0.00016663, 0.00033833, 0.00082988, 0.0034005, 0.0042941, 0.0050884, 0.0052612, 0.0055901, 0.0054357, 0.0052671, 0.0049174, 0.0042236, 0.0031138, 0.0011733, 0.0014057, 0.0010843, 0.0010992, 0.0007966, 0.00052196, 0.00053029, 0.00021273, 0.00041664, 0.00010455, 0.00015173, -9.7827e-05, -0.0010859, -0.0013748, -0.0016641, -0.0016887};
50 pt_density = p0[ind] + p1[ind]*(NPV - 1) + p2[ind]*(NPV - 1)*(NPV - 1);
52 double ECorr = pt_density*area*cosh(eta);
62 boost::regex re(
std::string(
"^(")+name+
"|"+name+
"_v\\d*)$");
63 for (
unsigned int i = 0,
n = list.size() ;
i <
n ; ++
i) {
64 if(boost::regex_match(list[
i],re))
return i;
91 if (iConfig.
exists(
"prodProcess"))
114 if (workOnAOD_ < 2) {
121 tok_HF_ = consumes<edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit> > >(
hfRecHitName_);
122 tok_HO_ = consumes<edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit> > >(
hoRecHitName_);
138 const char*
prod=
"GammaJetProd";
139 if (prodProcess_.size()==0) {
140 edm::LogError(
"GammaJetAnalysis") <<
"prodProcess needs to be defined";
143 const char* an= prodProcess_.c_str();
144 edm::LogWarning(
"GammaJetAnalysis") <<
"FAST FIX: changing " << photonCollName_
151 tok_HBHE_ = consumes<edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit> > >(
edm::InputTag(prod,hbheRecHitName_,an));
152 tok_HF_ = consumes<edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit> > >(
edm::InputTag(prod,hfRecHitName_,an));
153 tok_HO_ = consumes<edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit> > >(
edm::InputTag(prod,hoRecHitName_,an));
167 TString HLTlabel =
"TriggerResults::HLT";
168 if (prodProcess_.find(
"reRECO")!=std::string::npos)
169 HLTlabel.ReplaceAll(
"HLT",
"reHLT");
190 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find PhotonCollection named "
210 edm::LogWarning(
"GammaJetAnalysis") <<
"Failed to get photon quality flags";
217 if (!loosePhotonQual_Vec.
isValid() || !tightPhotonQual_Vec.
isValid()) {
218 edm::LogWarning(
"GammaJetAnalysis") <<
"Failed to get photon quality flags (vec)";
225 std::set<PhotonPair, PhotonPairComp> photonpairset;
227 for(reco::PhotonCollection::const_iterator it=photons->begin(); it!=photons->end(); ++it) {
236 HERE(Form(
"photonpairset.size=%d",
int(photonpairset.size())));
239 if (
debug_>0)
edm::LogInfo(
"GammaJetAnalysis") <<
"No good quality photons in the event";
251 for(std::set<PhotonPair, PhotonPairComp>::const_iterator it=photonpairset.begin(); it!=photonpairset.end(); ++it) {
254 if(counter==1) photon_tag = photon;
255 else if (counter==2) photon_2nd = photon;
264 HERE(Form(
"counter=%d",counter));
279 unsigned int anyJetCount=0;
287 anyJetCount+= pfjets->size();
291 HERE(Form(
"anyJetCount=%d",anyJetCount));
293 if (anyJetCount==0) {
320 if (!photonTrigFlag || !jetTrigFlag) {
324 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find TriggerResults::HLT";
331 const std::vector<std::string> *trNames= & evTrigNames.
triggerNames();
332 for (
size_t i=0;
i<trNames->size(); ++
i) {
333 if (trNames->at(
i).find(
"_Photon")!=std::string::npos) {
345 if (
id==evTrigNames.
size()) {
350 int fired= triggerResults->accept(
id);
351 if (fired) photonTrigFlag=
true;
363 if (
id==evTrigNames.
size()) {
368 int fired= triggerResults->accept(
id);
369 if (fired) jetTrigFlag=
true;
376 if (!photonTrigFlag && !jetTrigFlag) {
381 HERE(
"start isolation");
411 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find GenJet vector named "
427 if(!genEventInfoProduct.
isValid()){
433 if (genEventInfoProduct->hasBinningValues()) {
434 eventPtHat_ = genEventInfoProduct->binningValues()[0];
492 HERE(
"calc charged pfiso");
495 HERE(
"got isolation");
502 HERE(Form(
"workOnAOD_<2. loose photon qual size=%d",
int(loosePhotonQual->size())));
505 HERE(Form(
"got photon ref, photon_tag.idx()=%d",photon_tag.
idx()));
523 HERE(
"reset pho gen");
532 for (std::vector<reco::GenParticle>::const_iterator itmc=genparticles->begin();
533 itmc!=genparticles->end(); itmc++) {
534 if (itmc->status() == 1 && itmc->pdgId()==22) {
536 itmc->eta(),itmc->phi());
549 HERE(
"run over PFJets");
559 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find HBHERecHit named "
568 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find HFRecHit named "
577 edm::LogWarning(
"GammaJetAnalysis") <<
"Could not find HORecHit named "
582 HERE(
"get geometry");
600 edm::LogInfo(
"GammaJetAnalysis") << (*ith).id().ieta() <<
" "
601 << (*ith).id().iphi();
622 HERE(
"Get primary vertices");
633 for(std::vector<reco::Vertex>::const_iterator it=pv->begin(); it!=pv->end(); ++it){
634 if(!it->isFake() && it->ndof() > 4) ++
pf_NPV_;
637 HERE(
"get corrector");
643 std::set<PFJetCorretPair, PFJetCorretPairComp> pfjetcorretpairset;
644 for(reco::PFJetCollection::const_iterator it=pfjets->begin(); it!=pfjets->end(); ++it) {
647 if (
deltaR(photon_tag,jet)<0.5)
continue;
656 for(std::set<PFJetCorretPair, PFJetCorretPairComp>::const_iterator it=pfjetcorretpairset.begin(); it!=pfjetcorretpairset.end(); ++it) {
659 if(jet_cntr==1) pfjet_probe =
jet;
660 else if(jet_cntr==2) pf_2ndjet =
jet;
661 else if(jet_cntr==3) pf_3rdjet =
jet;
665 HERE(
"reached selection");
670 if (!pfjet_probe.
isValid()) failSelPF |= 1;
674 if (
deltaR(photon_tag,pfjet_probe.
jet())<0.5) failSelPF |= 4;
714 for (
int iJet=2; iJet>0; iJet--) {
718 if (iJet==2) pfjet_probe= pf_2ndjet;
719 else pfjet_probe = pfjet_probe_store;
721 if(!pfjet_probe.
jet()) {
726 edm::LogWarning(
"GammaJetAnalysis") <<
"error in the code: leading pf jet is null";
731 HERE(
"work further");
734 std::map<int,std::pair<int,std::set<float>>> ppfjet_rechits;
735 std::map<float,int> ppfjet_clusters;
756 HERE(
"Get PF constituents");
760 HERE(Form(
"probeconst.size=%d",
int(probeconst.size())));
762 for(std::vector<reco::PFCandidatePtr>::const_iterator it=probeconst.begin(); it!=probeconst.end(); ++it){
763 bool hasTrack =
false;
764 if (!(*it))
HERE(
"\tnull probeconst iterator value\n");
767 HERE(Form(
"iPF=%d",iPF));
770 switch(candidateType){
795 for(std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin(); itmc != genparticles->end(); itmc++){
796 if(itmc->status() == 1 && itmc->pdgId() > 100){
797 double dr =
deltaR((*it)->eta(),(*it)->phi(),itmc->eta(),itmc->phi());
800 genE = itmc->energy();
801 genpdgId = itmc->pdgId();
866 for(std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin(); itmc != genparticles->end(); itmc++){
867 if(itmc->status() == 1 && itmc->pdgId() > 100){
868 double dr =
deltaR((*it)->eta(),(*it)->phi(),itmc->eta(),itmc->phi());
871 genE = itmc->energy();
872 genpdgId = itmc->pdgId();
899 for(std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin(); itmc != genparticles->end(); itmc++){
900 if(itmc->status() == 1 && itmc->pdgId() > 100){
901 double dr =
deltaR((*it)->eta(),(*it)->phi(),itmc->eta(),itmc->phi());
904 genE = itmc->energy();
905 genpdgId = itmc->pdgId();
932 for(std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin(); itmc != genparticles->end(); itmc++){
933 if(itmc->status() == 1 && itmc->pdgId() > 100){
934 double dr =
deltaR((*it)->eta(),(*it)->phi(),itmc->eta(),itmc->phi());
937 genE = itmc->energy();
938 genpdgId = itmc->pdgId();
955 int maxElement=(*it)->elementsInBlocks().size();
961 HERE(Form(
"maxElement=%d",maxElement));
962 for(
int e=0;
e<maxElement; ++
e){
966 for(
unsigned iEle=0; iEle<elements.
size(); iEle++) {
967 if(elements[iEle].
index() == (*it)->elementsInBlocks()[
e].second){
974 if(ppfjet_clusters.count(cluster_dR) == 0){
981 int cluster_ind = ppfjet_clusters[cluster_dR];
982 std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.
hitsAndFractions();
985 int nHits = hitsAndFracs.size();
986 for(
int iHit=0; iHit<nHits; iHit++){
990 int etaPhiRecHit =
getEtaPhi((*ith).id());
991 if(etaPhiPF == etaPhiRecHit){
993 if(ppfjet_rechits.count((*ith).id()) == 0){
1000 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1010 switch((*ith).id().subdet()){
1015 float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1016 float avgphi = (
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1017 if((cv[0].
phi() < cv[2].phi()) && (
debug_>1))
edm::LogInfo(
"GammaJetAnalysis") <<
"pHB" << cv[0].phi() <<
" " << cv[2].phi();
1018 if(cv[0].
phi() < cv[2].
phi()) avgphi = (2.0*3.141592653 +
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1026 float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1027 float avgphi = (
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1028 if((cv[0].
phi() < cv[2].phi()) && (
debug_>1))
edm::LogInfo(
"GammaJetAnalysis") <<
"pHE" << cv[0].phi() <<
" " << cv[2].phi();
1029 if(cv[0].
phi() < cv[2].
phi()) avgphi = (2.0*3.141592653 +
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1040 else if(ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].
second) == 0){
1041 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1045 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1060 if((*ith).id().depth() == 1)
continue;
1064 bool passMatch =
false;
1065 if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
1066 if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch =
true;
1067 else if(cv[0].
phi() < cv[2].phi()){
1068 if((*it)->phi() < cv[0].phi()) passMatch =
true;
1069 else if((*it)->phi() > cv[2].phi()) passMatch =
true;
1085 float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1086 float avgphi = (
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1087 if((cv[0].
phi() < cv[2].phi()) && (
debug_>1))
edm::LogInfo(
"GammaJetAnalysis") <<
"pHFhad" << cv[0].phi() <<
" " << cv[2].phi();
1088 if(cv[0].
phi() < cv[2].
phi()) avgphi = (2.0*3.141592653 +
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1091 HFHAD_E += (*ith).energy();
1102 if((*ith).id().depth() == 2)
continue;
1106 bool passMatch =
false;
1107 if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
1108 if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch =
true;
1109 else if(cv[0].
phi() < cv[2].phi()){
1110 if((*it)->phi() < cv[0].phi()) passMatch =
true;
1111 else if((*it)->phi() > cv[2].phi()) passMatch =
true;
1127 float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1128 float avgphi = (
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1129 if((cv[0].
phi() < cv[2].phi()) && (
debug_>1))
edm::LogInfo(
"GammaJetAnalysis") <<
"pHFem" << cv[0].phi() <<
" " << cv[2].phi();
1130 if(cv[0].
phi() < cv[2].
phi()) avgphi = (2.0*3.141592653 +
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1133 HFEM_E += (*ith).energy();
1144 if(ppfjet_clusters.count(cluster_dR) == 0){
1151 int cluster_ind = ppfjet_clusters[cluster_dR];
1153 std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.
hitsAndFractions();
1154 int nHits = hitsAndFracs.size();
1155 for(
int iHit=0; iHit<nHits; iHit++){
1159 int etaPhiRecHit =
getEtaPhi((*ith).id());
1160 if(etaPhiPF == etaPhiRecHit){
1162 if(ppfjet_rechits.count((*ith).id()) == 0){
1169 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1181 float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1182 float avgphi = (
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1183 if((cv[0].
phi() < cv[2].phi()) && (
debug_>1))
edm::LogInfo(
"GammaJetAnalysis") <<
"pHO" << cv[0].phi() <<
" " << cv[2].phi();
1184 if(cv[0].
phi() < cv[2].
phi()) avgphi = (2.0*3.141592653 +
static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].
phi()))/2.0;
1190 else if(ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].
second) == 0){
1191 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1195 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1204 switch(candidateType){
1222 for(std::vector<reco::GenJet>::const_iterator it=genjets->begin(); it!=genjets->end(); ++it){
1224 double dr=
deltaR(jet, pfjet_probe.
jet());
1249 met_phi_ = pfmet_h->begin()->phi();
1254 if ( pfmetType1_h.
isValid()) {
1280 pf_tree_ =
new TTree(
"pf_gammajettree",
"tree for gamma+jet balancing using PFJets");
1284 for (
int iJet=0; iJet<2; iJet++) {
1286 if (!doJet)
continue;
1287 if (iJet>0)
continue;
1297 tree->Branch(
"RunNumber",&
runNumber_,
"RunNumber/I");
1298 tree->Branch(
"LumiBlock",&
lumiBlock_,
"LumiBlock/I");
1299 tree->Branch(
"EventNumber",&
eventNumber_,
"EventNumber/I");
1300 tree->Branch(
"EventWeight",&
eventWeight_,
"EventWeight/F");
1301 tree->Branch(
"EventPtHat",&
eventPtHat_,
"EventPtHat/F");
1304 tree->Branch(
"rho2012", &
rho2012_,
"rho2012/F");
1305 tree->Branch(
"tagPho_pt", &
tagPho_pt_,
"tagPho_pt/F");
1306 tree->Branch(
"pho_2nd_pt", &
pho_2nd_pt_,
"pho_2nd_pt/F");
1307 tree->Branch(
"tagPho_energy", &
tagPho_energy_,
"tagPho_energy/F");
1308 tree->Branch(
"tagPho_eta", &
tagPho_eta_,
"tagPho_eta/F");
1309 tree->Branch(
"tagPho_phi", &
tagPho_phi_,
"tagPho_phi/F");
1310 tree->Branch(
"tagPho_sieie", &
tagPho_sieie_,
"tagPho_sieie/F");
1311 tree->Branch(
"tagPho_HoE", &
tagPho_HoE_,
"tagPho_HoE/F");
1312 tree->Branch(
"tagPho_r9", &
tagPho_r9_,
"tagPho_r9/F");
1326 tree->Branch(
"tagPho_genPt",&
tagPho_genPt_,
"tagPho_genPt/F");
1328 tree->Branch(
"tagPho_genEta",&
tagPho_genEta_,
"tagPho_genEta/F");
1329 tree->Branch(
"tagPho_genPhi",&
tagPho_genPhi_,
"tagPho_genPhi/F");
1333 tree->Branch(
"nPhotons",&
nPhotons_,
"nPhotons/I");
1334 tree->Branch(
"nGenJets",&
nGenJets_,
"nGenJets/I");
1544 misc_tree_=
new TTree(
"misc_tree",
"tree for misc.info");
1555 TString str = TString(asctime(localtime(<ime)));
1556 if (str[str.Length()-1]==
'\n') str.Remove(str.Length()-1,1);
1557 TObjString date(str);
1558 date.Write(str.Data());
1577 if (!noPhotonTrigger &&
1580 if (!noJetTrigger &&
1583 if (noPhotonTrigger && noJetTrigger) {
1585 if (
debug_>1)
edm::LogInfo(
"GammaJetAnalysis") <<
"HLT trigger ignored: no trigger requested";
1594 if (
debug_>0)
edm::LogInfo(
"GammaJetAnalysis") <<
"Initializing trigger information for individual run";
1608 <<
" HLT config extraction failure with process name " <<
processName;
1624 if (localPho->
isEB()) {
1625 dRVeto = dRVetoBarrel;
1628 dRVeto = dRVetoEndcap;
1632 int nsize = forIsolation->size();
1634 for (
int i=0;
i<nsize;
i++) {
1643 if (localPho->
isEB()) {
1644 if (fabs(pfc.
pt()) < energyBarrel)
1647 if (fabs(pfc.
energy()) < energyEndcap)
1656 float dEta = fabs(photon_directionWrtVtx.Eta() - pfc.
momentum().Eta());
1657 float dR =
deltaR(photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.
momentum().Eta(), pfc.
momentum().Phi());
1659 if (dEta < etaStrip)
1662 if(dR > dRmax || dR < dRVeto)
1675 return pfEcalIso(localPho, pfHandle, dRmax, dRveto, dRveto, 0.0, 0.0, 0.0, 0.0, pfToUse);
1687 if (localPho->
isEB())
1688 dRveto = dRvetoBarrel;
1690 dRveto = dRvetoEndcap;
1692 std::vector<float>
result;
1696 if (
debug_>1)
edm::LogInfo(
"GammaJetAnalysis") <<
"vtxHandle->size() = " << vtxHandle->size();
1697 for(
unsigned int ivtx=0; ivtx<(vtxHandle->size()); ++ivtx) {
1704 if (
debug_>1)
edm::LogInfo(
"GammaJetAnalysis") <<
"pfTkIsoWithVertex :: Will Loop over the PFCandidates";
1707 for(
unsigned i=0;
i<forIsolation->size();
i++) {
1713 edm::LogInfo(
"GammaJetAnalysis") <<
"pfToUse=" << pfToUse;
1719 edm::LogInfo(
"GammaJetAnalysis") <<
"\n ***** HERE pfc.particleId() == pfToUse ";
1722 if (pfc.
pt() <
ptMin)
continue;
1724 float dz = fabs(pfc.
trackRef()->dz(vtx->position()));
1725 if (dz > dzMax)
continue;
1727 float dxy = fabs(pfc.
trackRef()->dxy(vtx->position()));
1728 if(fabs(dxy) > dxyMax)
continue;
1729 float dR =
deltaR(photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.
momentum().Eta(), pfc.
momentum().Phi());
1730 if(dR > dRmax || dR < dRveto)
continue;
1737 result.push_back(sum);
1740 edm::LogInfo(
"GammaJetAnalysis") <<
"Will return result";
1899 double deta = j1->
eta()-j2->
eta();
1900 double dphi = std::fabs(j1->
phi()-j2->
phi());
1901 if(dphi>3.1415927) dphi = 2*3.1415927 - dphi;
1902 return std::sqrt(deta*deta + dphi*dphi);
1909 double deta = eta1 - eta2;
1910 double dphi = std::fabs(phi1 - phi2);
1911 if(dphi>3.1415927) dphi = 2*3.1415927 - dphi;
1912 return std::sqrt(deta*deta + dphi*dphi);
1932 return id.rawId() & 0x3FFF;
1939 return id.rawId() & 0x3FFF;
float tagPho_pfiso_myphoton03_
T getParameter(std::string const &) const
EventNumber_t event() const
virtual Photon * clone() const
returns a clone of the candidate
T getUntrackedParameter(std::string const &, T const &) const
void copy_leadingPfJetVars_to_pfJet2()
std::vector< float > ppfjet_cluster_phi_
float hcalTowerSumEtConeDR04() const
Hcal isolation sum.
float tagPho_TrkIsoHollowDR04_
virtual double p() const
magnitude of momentum vector
std::vector< int > pfjet2_twr_ieta_
std::vector< int > pfjet2_had_id_
edm::EDGetTokenT< std::vector< Bool_t > > tok_tightPhotonV_
std::vector< int > pfjet2_twr_iphi_
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
float tagPho_HcalIsoDR04_
ParticleType
particle types
std::vector< float > pfjet2_candtrack_px_
bool isNonnull() const
Checks for non-null.
std::vector< float > ppfjet_had_py_
std::vector< float > ppfjet_candtrack_px_
edm::EDGetTokenT< reco::ConversionCollection > tok_Conv_
edm::EDGetTokenT< GenEventInfoProduct > tok_GenEvInfo_
int tagPho_ConvSafeEleVeto_
float ppfjet_electron_pz_
reco::SuperClusterRef superCluster() const
Ref to SuperCluster.
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
float pfjet2_electron_EcalE_
virtual double et() const
transverse energy
std::vector< int > pfjet2_twr_hadind_
std::vector< float > ppfjet_candtrack_pz_
std::vector< int > pfjet2_twr_elmttype_
edm::EDGetTokenT< reco::PFCandidateCollection > tok_PFCand_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::vector< int > photonTrigFired_
std::vector< float > pfjet2_cluster_dR_
virtual double correction(const LorentzVector &fJet) const =0
get correction using Jet information only
#define DEFINE_FWK_MODULE(type)
edm::EDGetTokenT< edm::SortedCollection< HFRecHit, edm::StrictWeakOrdering< HFRecHit > > > tok_HF_
float chargedHadronEnergyFraction() const
chargedHadronEnergyFraction
Base class for all types of Jets.
float ppfjet_ChargedEMFrac_
std::vector< float > ppfjet_had_E_
float pfjet2_NeutralEMFrac_
edm::EDGetTokenT< double > tok_Rho_
edm::EDGetTokenT< edm::TriggerResults > tok_TrigRes_
std::vector< int > ppfjet_twr_ieta_
std::vector< int > ppfjet_twr_hadind_
std::string hoRecHitName_
bool writeTriggerPrescale_
std::vector< int > ppfjet_twr_subdet_
virtual Vector momentum() const
spatial momentum vector
int ppfjet_nConstituents_
bool exists(std::string const ¶meterName) const
checks if a parameter exists
std::vector< int > ppfjet_had_candtrackind_
int pfjet2_nConstituents_
float ecalRecHitSumEtConeDR04() const
const reco::PFJet * jet(void) const
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
std::vector< float > pfTkIsoWithVertex(const reco::Photon *localPho1, edm::Handle< reco::PFCandidateCollection > pfHandle, edm::Handle< reco::VertexCollection > vtxHandle, float dRmax, float dRvetoBarrel, float dRvetoEndcap, float ptMin, float dzMax, float dxyMax, reco::PFCandidate::ParticleType pfToUse)
std::string hbheRecHitName_
std::vector< float > pfjet2_had_E_
Strings::size_type size() const
float ppfjet_NeutralHadronFrac_
edm::EDGetTokenT< reco::PFMETCollection > tok_PFMET_
edm::EDGetTokenT< reco::VertexCollection > tok_Vertex_
double px() const
x coordinate of momentum vector
float ppfjet_electron_EcalE_
std::vector< float > pfjet2_had_EcalE_
float ppfjet_ChargedHadronFrac_
std::vector< float > pfjet2_cluster_eta_
LuminosityBlockNumber_t luminosityBlock() const
float tagPho_EcalIsoDR04_
edm::EDGetTokenT< reco::BeamSpot > tok_BS_
int chargedMultiplicity() const
chargedMultiplicity
std::vector< float > pfjet2_candtrack_py_
Strings const & triggerNames() const
Jets made from PFObjects.
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
std::vector< int > pfjet2_twr_candtrackind_
double eta() const
pseudorapidity of cluster centroid
std::string photonCollName_
std::vector< int > pfjet2_twr_clusterind_
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< int > pfjet2_twr_depth_
std::vector< int > pfjet2_had_candtrackind_
std::vector< float > ppfjet_twr_frac_
reco::TrackRef trackRef() const
U second(std::pair< T, U > const &p)
float pfjet2_electron_px_
edm::EDGetTokenT< std::vector< reco::GenParticle > > tok_GenPart_
float pfjet2_NeutralHadronFrac_
std::vector< int > ppfjet_had_id_
float pfjet2_ChargedHadronFrac_
float neutralEmEnergyFraction() const
neutralEmEnergyFraction
virtual double energy() const
energy
float met_value_
MET info.
std::vector< int > jetTrigPrescale_
std::vector< int > ppfjet_twr_clusterind_
edm::EDGetTokenT< edm::SortedCollection< HORecHit, edm::StrictWeakOrdering< HORecHit > > > tok_HO_
float ppfjet_NeutralEMFrac_
std::string rootHistFilename_
float pfjet2_electron_pz_
std::string hfRecHitName_
std::vector< float > pfjet2_had_E_mctruth_
std::vector< float > ppfjet_candtrack_py_
float pfEcalIso(const reco::Photon *localPho1, edm::Handle< reco::PFCandidateCollection > pfHandle, float dRmax, float dRVetoBarrel, float dRVetoEndcap, float etaStripBarrel, float etaStripEndcap, float energyBarrel, float energyEndcap, reco::PFCandidate::ParticleType pfToUse)
std::vector< float > pfjet2_had_py_
float pfHcalIso(const reco::Photon *localPho, edm::Handle< reco::PFCandidateCollection > pfHandle, float dRmax, float dRveto, reco::PFCandidate::ParticleType pfToUse)
std::vector< float > ppfjet_cluster_dR_
edm::EDGetTokenT< reco::PhotonCollection > tok_Photon_
float pfjet2_ChargedEMFrac_
float ppfjet_ChargedMultiplicity_
float tagPho_pfiso_myneutral03_
float sigmaIetaIeta() const
std::vector< float > ppfjet_twr_hade_
std::vector< int > ppfjet_twr_elmttype_
unsigned int helper_findTrigger(const std::vector< std::string > &list, const std::string &name)
std::vector< float > ppfjet_had_rawHcalE_
edm::InputTag rhoCollection_
edm::EDGetTokenT< std::vector< Bool_t > > tok_loosePhotonV_
std::vector< int > ppfjet_twr_depth_
edm::EDGetTokenT< edm::ValueMap< Bool_t > > tok_tightPhoton_
float pfjet2_photon_EcalE_
const reco::Photon * photon(void) const
virtual const Point & vertex() const
vertex position (overwritten by PF...)
Jets made from MC generator particles.
float neutralHadronEnergyFraction() const
neutralHadronEnergyFraction
std::vector< float > pfjet2_candtrack_pz_
std::vector< float > pfjet2_twr_frac_
float hadTowOverEm() const
the ration of hadronic energy in towers behind the BCs in the SC and the SC energy ...
edm::InputTag pfType1METColl
std::vector< float > pfjet2_had_rawHcalE_
HLTConfigProvider hltConfig_
static std::string const triggerResults
float hadronicOverEm() const
the total hadronic over electromagnetic fraction
void beginRun(const edm::Run &, const edm::EventSetup &)
std::string genParticleCollName_
std::vector< float > ppfjet_twr_dR_
std::vector< float > ppfjet_candtrack_EcalE_
std::vector< float > ppfjet_cluster_eta_
double pz() const
z coordinate of momentum vector
float pfjet2_unkown_EcalE_
std::vector< float > pfjet2_cluster_phi_
std::string pfJetCollName_
double getNeutralPVCorr(double eta, int intNPV, double area, bool isMC_)
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
void clear_leadingPfJetVars()
float chargedEmEnergyFraction() const
chargedEmEnergyFraction
std::vector< std::string > photonTrigNamesV_
float ppfjet_photon_EcalE_
T const * product() const
virtual double px() const
x coordinate of momentum vector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
std::string const & triggerName(unsigned int index) const
std::vector< int > jetTrigFired_
XYZPointD XYZPoint
point in space with cartesian internal representation
edm::EDGetTokenT< reco::PFJetCollection > tok_PFJet_
std::vector< float > ppfjet_had_E_mctruth_
edm::EDGetTokenT< edm::ValueMap< Bool_t > > tok_loosePhoton_
bool hasPixelSeed() const
Bool flagging photons having a non-zero size vector of Ref to electornPixel seeds.
std::vector< std::string > jetTrigNamesV_
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d'tor
static const JetCorrector * getJetCorrector(const std::string &fName, const edm::EventSetup &fSetup)
retrieve corrector from the event setup. troughs exception if something is missing ...
edm::EDGetTokenT< std::vector< reco::Vertex > > tok_PV_
std::pair< int, int > prescaleValues(const edm::Event &iEvent, const edm::EventSetup &iSetup, const std::string &trigger) const
Combined L1T (pair.first) and HLT (pair.second) prescales per HLT path.
float ppfjet_electron_py_
std::vector< int > pfjet2_had_ntwrs_
std::string genEventInfoName_
double deltaR(const reco::Jet *j1, const reco::Jet *j2)
int getEtaPhi(const DetId id)
float trkSumPtHollowConeDR04() const
std::vector< float > pfjet2_twr_hade_
edm::EDGetTokenT< edm::SortedCollection< HBHERecHit, edm::StrictWeakOrdering< HBHERecHit > > > tok_HBHE_
std::string genJetCollName_
void HERE(const char *msg)
Particle reconstructed by the particle flow algorithm.
std::vector< float > pfjet2_had_emf_
virtual float jetArea() const
get jet area
virtual int nConstituents() const
of constituents
std::vector< float > ppfjet_had_px_
std::vector< int > ppfjet_had_mcpdgId_
static std::atomic< unsigned int > counter
std::vector< float > pfjet2_twr_dR_
GammaJetAnalysis(const edm::ParameterSet &)
float pfjet2_electron_py_
std::vector< float > pfjet2_had_pz_
edm::EDGetTokenT< reco::GsfElectronCollection > tok_GsfElec_
float pfjet2_ChargedMultiplicity_
std::vector< int > pfjet2_twr_subdet_
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
std::vector< int > ppfjet_twr_iphi_
std::vector< int > photonTrigPrescale_
std::vector< int > ppfjet_had_ntwrs_
edm::EDGetTokenT< std::vector< reco::GenJet > > tok_GenJet_
std::vector< float > ppfjet_had_pz_
std::vector< float > pfjet2_candtrack_EcalE_
float ppfjet_unkown_EcalE_
std::vector< int > pfjet2_had_mcpdgId_
edm::EDGetTokenT< reco::PFMETCollection > tok_PFType1MET_
float ppfjet_electron_px_
double phi() const
azimuthal angle of cluster centroid
virtual ParticleType particleId() const
virtual void analyze(const edm::Event &, const edm::EventSetup &)
const CornersVec & getCorners() const
Returns the corner points of this cell's volume.
std::vector< float > ppfjet_had_emf_
virtual double phi() const
momentum azimuthal angle
std::vector< int > ppfjet_twr_candtrackind_
float calc_dPhi(const PhotonPair &pho, const JetPair_type &jet)
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
virtual double py() const
y coordinate of momentum vector
std::vector< std::vector< float > > tagPho_pfiso_mycharged03
std::vector< float > ppfjet_had_EcalE_
std::vector< float > pfjet2_had_px_
double py() const
y coordinate of momentum vector
std::string pfJetCorrName_
float tagPho_HcalIsoDR0412_
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger list("!*","!HLTx*"if it matches 2 triggers or more) will accept the event if all the matching triggers are FAIL.It will reject the event if any of the triggers are PASS or EXCEPTION(this matches the behavior of"!*"before the partial wildcard feature was incorporated).Triggers which are in the READY state are completely ignored.(READY should never be returned since the trigger paths have been run
reco::SuperClusterRef superClusterRef() const
return a reference to the corresponding SuperCluster if any