| 1 |
arizzi |
1.1 |
#include <TH1F.h>
|
| 2 |
|
|
#include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
|
| 3 |
|
|
#include <TH2F.h>
|
| 4 |
|
|
#include <TROOT.h>
|
| 5 |
|
|
#include <TFile.h>
|
| 6 |
|
|
#include <TTree.h>
|
| 7 |
|
|
#include <TSystem.h>
|
| 8 |
|
|
#include "DataFormats/FWLite/interface/Event.h"
|
| 9 |
|
|
#include "DataFormats/FWLite/interface/Handle.h"
|
| 10 |
|
|
#include "FWCore/FWLite/interface/AutoLibraryLoader.h"
|
| 11 |
|
|
#include "DataFormats/MuonReco/interface/Muon.h"
|
| 12 |
|
|
#include "DataFormats/PatCandidates/interface/Muon.h"
|
| 13 |
|
|
#include "PhysicsTools/FWLite/interface/TFileService.h"
|
| 14 |
|
|
#include "FWCore/ParameterSet/interface/ProcessDesc.h"
|
| 15 |
|
|
#include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
|
| 16 |
|
|
#include "DataFormats/Math/interface/deltaR.h"
|
| 17 |
|
|
#include "DataFormats/Math/interface/deltaPhi.h"
|
| 18 |
|
|
|
| 19 |
|
|
#include "DataFormats/FWLite/interface/LuminosityBlock.h"
|
| 20 |
|
|
#include "DataFormats/FWLite/interface/Run.h"
|
| 21 |
|
|
#include "DataFormats/Luminosity/interface/LumiSummary.h"
|
| 22 |
|
|
|
| 23 |
arizzi |
1.4 |
#include "VHbbAnalysis/VHbbDataFormats/interface/HbbCandidateFinderAlgo.h"
|
| 24 |
|
|
#include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc"
|
| 25 |
arizzi |
1.1 |
|
| 26 |
|
|
#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
|
| 27 |
|
|
#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
|
| 28 |
|
|
#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbCandidate.h"
|
| 29 |
|
|
#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerReader.h"
|
| 30 |
nmohr |
1.5 |
#include "VHbbAnalysis/VHbbDataFormats/interface/TopMassReco.h"
|
| 31 |
arizzi |
1.1 |
|
| 32 |
bortigno |
1.30 |
//for IVF
|
| 33 |
|
|
#include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
|
| 34 |
|
|
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
|
| 35 |
|
|
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
|
| 36 |
|
|
#include "DataFormats/GeometryVector/interface/VectorUtil.h"
|
| 37 |
|
|
#include <DataFormats/GeometrySurface/interface/Surface.h>
|
| 38 |
|
|
#include "Math/SMatrix.h"
|
| 39 |
arizzi |
1.16 |
|
| 40 |
|
|
//Move class definition to Ntupler.h ?
|
| 41 |
|
|
//#include "VHbbAnalysis/VHbbDataFormats/interface/Ntupler.h"
|
| 42 |
|
|
|
| 43 |
bortigno |
1.30 |
//for SimBhad: FIXME
|
| 44 |
|
|
//#include "VHbbAnalysis/BAnalysis/interface/SimBHadron.h"
|
| 45 |
|
|
//btagging
|
| 46 |
arizzi |
1.16 |
#include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
|
| 47 |
bortigno |
1.30 |
//trigger
|
| 48 |
arizzi |
1.16 |
#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
|
| 49 |
|
|
|
| 50 |
arizzi |
1.1 |
#include <sstream>
|
| 51 |
|
|
#include <string>
|
| 52 |
|
|
|
| 53 |
|
|
#define MAXJ 30
|
| 54 |
|
|
#define MAXL 10
|
| 55 |
bortigno |
1.30 |
#define MAXB 10
|
| 56 |
arizzi |
1.1 |
|
| 57 |
arizzi |
1.24 |
struct CompareDeltaR {
|
| 58 |
bortigno |
1.30 |
CompareDeltaR(TLorentzVector p4dir_): p4dir(p4dir_) {}
|
| 59 |
|
|
bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
|
| 60 |
|
|
return j1.p4.DeltaR(p4dir) > j2.p4.DeltaR(p4dir);
|
| 61 |
|
|
}
|
| 62 |
|
|
TLorentzVector p4dir;
|
| 63 |
arizzi |
1.24 |
};
|
| 64 |
arizzi |
1.4 |
|
| 65 |
bortigno |
1.30 |
const GlobalVector flightDirection(const TVector3 pv, const reco::Vertex &sv){
|
| 66 |
|
|
GlobalVector fdir(sv.position().X() - pv.X(),
|
| 67 |
|
|
sv.position().Y() - pv.Y(),
|
| 68 |
|
|
sv.position().Z() - pv.Z());
|
| 69 |
|
|
return fdir;
|
| 70 |
|
|
}
|
| 71 |
|
|
|
| 72 |
nmohr |
1.7 |
bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
|
| 73 |
|
|
const edm::EventBase &event)
|
| 74 |
|
|
{
|
| 75 |
bortigno |
1.30 |
// if the jsonVec is empty, then no JSON file was provided so all
|
| 76 |
|
|
// events should pass
|
| 77 |
|
|
if (jsonVec.empty())
|
| 78 |
|
|
{
|
| 79 |
nmohr |
1.7 |
return true;
|
| 80 |
bortigno |
1.30 |
}
|
| 81 |
|
|
bool (* funcPtr) (edm::LuminosityBlockRange const &,
|
| 82 |
|
|
edm::LuminosityBlockID const &) = &edm::contains;
|
| 83 |
|
|
edm::LuminosityBlockID lumiID (event.id().run(),
|
| 84 |
|
|
event.id().luminosityBlock());
|
| 85 |
|
|
std::vector< edm::LuminosityBlockRange >::const_iterator iter =
|
| 86 |
|
|
std::find_if (jsonVec.begin(), jsonVec.end(),
|
| 87 |
|
|
boost::bind(funcPtr, _1, lumiID) );
|
| 88 |
|
|
return jsonVec.end() != iter;
|
| 89 |
nmohr |
1.7 |
|
| 90 |
|
|
}
|
| 91 |
|
|
|
| 92 |
arizzi |
1.1 |
|
| 93 |
bortigno |
1.30 |
//FIXME : need to update EDM ntuple with SimBhadron infos
|
| 94 |
|
|
// typedef struct
|
| 95 |
|
|
// {
|
| 96 |
|
|
// void set(const SimBHadron & sbhc, int i){
|
| 97 |
|
|
// mass[i] = sbhc.mass();
|
| 98 |
|
|
// pt[i] = sbhc.pt();
|
| 99 |
|
|
// eta[i] = sbhc.eta();
|
| 100 |
|
|
// phi[i] = sbhc.phi();
|
| 101 |
|
|
// vtx_x[i] = sbhc.decPosition.x();
|
| 102 |
|
|
// vtx_y[i] = sbhc.decPosition.y();
|
| 103 |
|
|
// vtx_z[i] = sbhc.decPosition.z();
|
| 104 |
|
|
// pdgId[i] = sbhc.pdgId();
|
| 105 |
|
|
// status[i] = sbhc.status();
|
| 106 |
|
|
// };
|
| 107 |
|
|
// void reset(){
|
| 108 |
|
|
// for(int i=0; i < MAXB; ++i){
|
| 109 |
|
|
// mass[i] = -99; pt[i] = -99; eta[i] = -99; phi[i] = -99; vtx_x[i] = -99; vtx_y[i] = -99; vtx_z[i] = -99; pdgId[i] = -99; status[i] = -99;
|
| 110 |
|
|
// }
|
| 111 |
|
|
// };
|
| 112 |
|
|
// float mass[MAXB];
|
| 113 |
|
|
// float pt[MAXB];
|
| 114 |
|
|
// float eta[MAXB];
|
| 115 |
|
|
// float phi[MAXB];
|
| 116 |
|
|
// float vtx_x[MAXB];
|
| 117 |
|
|
// float vtx_y[MAXB];
|
| 118 |
|
|
// float vtx_z[MAXB];
|
| 119 |
|
|
// int pdgId[MAXB];
|
| 120 |
|
|
// int status[MAXB];
|
| 121 |
|
|
// // int quarkStatus[MAXB];
|
| 122 |
|
|
// // int brotherStatus[MAXB];
|
| 123 |
|
|
// // int otherId[MAXB];
|
| 124 |
|
|
// // bool etaOk[MAXB];
|
| 125 |
|
|
// // bool simOk[MAXB];
|
| 126 |
|
|
// // bool trackOk[MAXB];
|
| 127 |
|
|
// // bool cutOk[MAXB];
|
| 128 |
|
|
// // bool cutNewOk[MAXB];
|
| 129 |
|
|
// // bool mcMatchOk[MAXB];
|
| 130 |
|
|
// // bool matchOk[MAXB];
|
| 131 |
|
|
// } SimBHadronInfo;
|
| 132 |
|
|
|
| 133 |
|
|
|
| 134 |
|
|
typedef struct
|
| 135 |
|
|
{
|
| 136 |
|
|
void set( const reco::SecondaryVertex & recoSv, const TVector3 recoPv, int isv){
|
| 137 |
|
|
pt[isv] = recoSv.p4().Pt();
|
| 138 |
|
|
eta[isv] = flightDirection(recoPv,recoSv).eta();
|
| 139 |
|
|
phi[isv] = flightDirection(recoPv,recoSv).phi();
|
| 140 |
|
|
massBcand[isv] = recoSv.p4().M();
|
| 141 |
|
|
massSv[isv] = recoSv.p4().M();
|
| 142 |
|
|
dist3D[isv] = recoSv.dist3d().value();
|
| 143 |
|
|
distSig3D[isv] = recoSv.dist3d().significance();
|
| 144 |
|
|
dist2D[isv] = recoSv.dist2d().value();
|
| 145 |
|
|
distSig2D[isv] = recoSv.dist2d().significance();
|
| 146 |
|
|
dist3D_norm[isv] = recoSv.dist3d().value()/recoSv.p4().Gamma();
|
| 147 |
|
|
};
|
| 148 |
|
|
void reset(){
|
| 149 |
|
|
for(int i = 0; i < MAXB; ++i){
|
| 150 |
|
|
massBcand[i] = -99; massSv[i]= -99; pt[i] = -99; eta[i] = -99; phi[i] = -99; dist3D[i] = -99; distSig3D[i] = -99; dist2D[i] = -99; distSig2D[i] = -99; dist3D_norm[i] = -99;
|
| 151 |
|
|
}
|
| 152 |
|
|
};
|
| 153 |
|
|
float massBcand[MAXB];
|
| 154 |
|
|
float massSv[MAXB];
|
| 155 |
|
|
float pt[MAXB];
|
| 156 |
|
|
float eta[MAXB];
|
| 157 |
|
|
float phi[MAXB];
|
| 158 |
|
|
float dist3D[MAXB];
|
| 159 |
|
|
float distSig3D[MAXB];
|
| 160 |
|
|
float dist2D[MAXB];
|
| 161 |
|
|
float distSig2D[MAXB];
|
| 162 |
|
|
float dist3D_norm[MAXB];
|
| 163 |
|
|
} IVFInfo;
|
| 164 |
arizzi |
1.1 |
|
| 165 |
bortigno |
1.30 |
|
| 166 |
|
|
typedef struct
|
| 167 |
|
|
{
|
| 168 |
|
|
float mass; //MT in case of W
|
| 169 |
|
|
float pt;
|
| 170 |
|
|
float eta;
|
| 171 |
|
|
float phi;
|
| 172 |
|
|
} TrackInfo;
|
| 173 |
arizzi |
1.1 |
|
| 174 |
|
|
|
| 175 |
arizzi |
1.9 |
struct LeptonInfo
|
| 176 |
bortigno |
1.30 |
{
|
| 177 |
|
|
void reset()
|
| 178 |
arizzi |
1.1 |
{
|
| 179 |
|
|
for(int i =0; i < MAXL;i++){ mass[i]=-99; pt[i]=-99; eta[i]=-99; phi[i]=-99; aodCombRelIso[i]=-99; pfCombRelIso[i]=-99; photonIso[i]=-99; neutralHadIso[i]=-99; chargedHadIso[i]=-99; particleIso[i]=-99; dxy[i]=-99; dz[i]=-99; type[i]=-99; }
|
| 180 |
bortigno |
1.30 |
}
|
| 181 |
arizzi |
1.1 |
|
| 182 |
bortigno |
1.30 |
template <class Input> void set(const Input & i, int j,int t)
|
| 183 |
|
|
{
|
| 184 |
|
|
type[j]=t;
|
| 185 |
|
|
pt[j]=i.p4.Pt();
|
| 186 |
|
|
mass[j]=i.p4.M();
|
| 187 |
|
|
eta[j]=i.p4.Eta();
|
| 188 |
|
|
phi[j]=i.p4.Phi();
|
| 189 |
|
|
aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
|
| 190 |
|
|
pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
|
| 191 |
|
|
photonIso[j]=i.pfPhoIso;
|
| 192 |
|
|
neutralHadIso[j]=i.pfNeuIso;
|
| 193 |
|
|
chargedHadIso[j]=i.pfChaIso;
|
| 194 |
|
|
setSpecific(i,j);
|
| 195 |
|
|
}
|
| 196 |
|
|
template <class Input> void setSpecific(const Input & i, int j)
|
| 197 |
|
|
{
|
| 198 |
|
|
}
|
| 199 |
arizzi |
1.4 |
|
| 200 |
|
|
|
| 201 |
|
|
|
| 202 |
arizzi |
1.1 |
|
| 203 |
bortigno |
1.30 |
float mass[MAXL]; //MT in case of W
|
| 204 |
|
|
float pt[MAXL];
|
| 205 |
|
|
float eta[MAXL];
|
| 206 |
|
|
float phi[MAXL];
|
| 207 |
|
|
float aodCombRelIso[MAXL];
|
| 208 |
|
|
float pfCombRelIso[MAXL];
|
| 209 |
|
|
float photonIso[MAXL];
|
| 210 |
|
|
float neutralHadIso[MAXL];
|
| 211 |
|
|
float chargedHadIso[MAXL];
|
| 212 |
|
|
float particleIso[MAXL];
|
| 213 |
|
|
float dxy[MAXL];
|
| 214 |
|
|
float dz[MAXL];
|
| 215 |
|
|
int type[MAXL];
|
| 216 |
|
|
float id80[MAXL];
|
| 217 |
|
|
float id95[MAXL];
|
| 218 |
|
|
};
|
| 219 |
arizzi |
1.1 |
|
| 220 |
bortigno |
1.30 |
template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
|
| 221 |
|
|
id80[j]=i.id80r;
|
| 222 |
|
|
id95[j]=i.id95r;
|
| 223 |
|
|
}
|
| 224 |
|
|
template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
|
| 225 |
|
|
dxy[j]=i.ipDb;
|
| 226 |
|
|
dz[j]=i.zPVPt;
|
| 227 |
|
|
}
|
| 228 |
|
|
|
| 229 |
arizzi |
1.4 |
|
| 230 |
bortigno |
1.30 |
typedef struct
|
| 231 |
|
|
{
|
| 232 |
|
|
float et;
|
| 233 |
|
|
float sumet;
|
| 234 |
|
|
float sig;
|
| 235 |
|
|
float phi;
|
| 236 |
|
|
} METInfo;
|
| 237 |
arizzi |
1.1 |
|
| 238 |
bortigno |
1.30 |
typedef struct
|
| 239 |
|
|
{
|
| 240 |
|
|
float mht;
|
| 241 |
|
|
float ht;
|
| 242 |
|
|
float sig;
|
| 243 |
|
|
float phi;
|
| 244 |
|
|
} MHTInfo;
|
| 245 |
nmohr |
1.5 |
|
| 246 |
bortigno |
1.30 |
typedef struct
|
| 247 |
|
|
{
|
| 248 |
|
|
float mass;
|
| 249 |
|
|
float pt;
|
| 250 |
|
|
float wMass;
|
| 251 |
|
|
} TopInfo;
|
| 252 |
arizzi |
1.1 |
|
| 253 |
bortigno |
1.30 |
typedef struct
|
| 254 |
|
|
{
|
| 255 |
|
|
int run;
|
| 256 |
|
|
int lumi;
|
| 257 |
|
|
int event;
|
| 258 |
|
|
int json;
|
| 259 |
|
|
} EventInfo;
|
| 260 |
arizzi |
1.1 |
|
| 261 |
bortigno |
1.30 |
typedef struct
|
| 262 |
|
|
{
|
| 263 |
|
|
void set(const VHbbEvent::SimpleJet & j, int i)
|
| 264 |
arizzi |
1.1 |
{
|
| 265 |
bortigno |
1.30 |
pt[i]=j.p4.Pt();
|
| 266 |
|
|
eta[i]=j.p4.Eta();
|
| 267 |
|
|
phi[i]=j.p4.Phi();
|
| 268 |
|
|
csv[i]=j.csv;
|
| 269 |
|
|
numTracksSV[i] = j.vtxMass;
|
| 270 |
|
|
vtxMass[i]= j.vtxMass;
|
| 271 |
|
|
vtx3dL[i] = j.vtx3dL;
|
| 272 |
|
|
vtx3deL[i] = j.vtx3deL;
|
| 273 |
|
|
chf[i]=j.chargedHadronEFraction;
|
| 274 |
|
|
nhf[i] =j.neutralHadronEFraction;
|
| 275 |
|
|
cef[i] =j.chargedEmEFraction;
|
| 276 |
|
|
nef[i] =j.neutralEmEFraction;
|
| 277 |
|
|
nconstituents[i] =j.nConstituents;
|
| 278 |
|
|
nch[i]=j.ntracks;
|
| 279 |
|
|
|
| 280 |
|
|
flavour[i]=j.flavour;
|
| 281 |
|
|
if(j.bestMCp4.Pt() > 0)
|
| 282 |
|
|
{
|
| 283 |
|
|
genPt[i]=j.bestMCp4.Pt();
|
| 284 |
|
|
genEta[i]=j.bestMCp4.Eta();
|
| 285 |
|
|
genPhi[i]=j.bestMCp4.Phi();
|
| 286 |
|
|
}
|
| 287 |
|
|
JECUnc[i]=j.jecunc;
|
| 288 |
arizzi |
1.3 |
|
| 289 |
bortigno |
1.30 |
}
|
| 290 |
|
|
void reset()
|
| 291 |
|
|
{
|
| 292 |
|
|
for(int i=0;i<MAXJ;i++) {
|
| 293 |
|
|
pt[i]=-99; eta[i]=-99; phi[i]=-99; csv[i]=-99; cosTheta[i]=-99; numTracksSV[i]=-99; chf[i]=-99; nhf[i]=-99; cef[i]=-99; nef[i]=-99; nch[i]=-99; nconstituents[i]=-99; flavour[i]=-99; genPt[i]=-99; genEta[i]=-99; genPhi[i]=-99; JECUnc[i]=-99;
|
| 294 |
arizzi |
1.1 |
}
|
| 295 |
bortigno |
1.30 |
}
|
| 296 |
|
|
float pt[MAXJ];
|
| 297 |
|
|
float eta[MAXJ];
|
| 298 |
|
|
float phi[MAXJ];
|
| 299 |
|
|
float csv[MAXJ];
|
| 300 |
|
|
float cosTheta[MAXJ];
|
| 301 |
|
|
int numTracksSV[MAXJ];
|
| 302 |
|
|
float chf[MAXJ];
|
| 303 |
|
|
float nhf[MAXJ];
|
| 304 |
|
|
float cef[MAXJ];
|
| 305 |
|
|
float nef[MAXJ];
|
| 306 |
|
|
float nch[MAXJ];
|
| 307 |
|
|
float nconstituents[MAXJ];
|
| 308 |
|
|
float flavour[MAXJ];
|
| 309 |
|
|
float genPt[MAXJ];
|
| 310 |
|
|
float genEta[MAXJ];
|
| 311 |
|
|
float genPhi[MAXJ];
|
| 312 |
|
|
float JECUnc[MAXJ];
|
| 313 |
|
|
float vtxMass[MAXJ];
|
| 314 |
|
|
float vtx3dL [MAXJ];
|
| 315 |
|
|
float vtx3deL[MAXJ];
|
| 316 |
|
|
} JetInfo;
|
| 317 |
arizzi |
1.1 |
|
| 318 |
|
|
int main(int argc, char* argv[])
|
| 319 |
|
|
{
|
| 320 |
|
|
gROOT->Reset();
|
| 321 |
|
|
|
| 322 |
|
|
TTree *_outTree;
|
| 323 |
bortigno |
1.30 |
IVFInfo IVF;
|
| 324 |
|
|
//FIXME
|
| 325 |
|
|
// SimBHadronInfo SimBs;
|
| 326 |
arizzi |
1.9 |
METInfo MET;
|
| 327 |
|
|
MHTInfo MHT;
|
| 328 |
|
|
TopInfo top;
|
| 329 |
|
|
EventInfo EVENT;
|
| 330 |
bortigno |
1.30 |
// JetInfo jet1,jet2, addJet1, addJet2;
|
| 331 |
arizzi |
1.9 |
// lepton1,lepton2;
|
| 332 |
|
|
JetInfo hJets, aJets;
|
| 333 |
|
|
LeptonInfo vLeptons, aLeptons;
|
| 334 |
arizzi |
1.1 |
int naJets=0, nhJets=0;
|
| 335 |
arizzi |
1.9 |
TrackInfo H;
|
| 336 |
|
|
TrackInfo V;
|
| 337 |
|
|
int nvlep=0,nalep=0;
|
| 338 |
arizzi |
1.1 |
|
| 339 |
bortigno |
1.30 |
float jjdr,jjdPhi,jjdEta,HVdPhi,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMET, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight;
|
| 340 |
|
|
float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart;
|
| 341 |
arizzi |
1.24 |
|
| 342 |
bortigno |
1.30 |
float BCands_dr = -99, BCands_dEta= -99, BCands_dPhi = -99, BCands_Hmass = -99, BCands_Hpt = -99, BCands_Heta = -99, BCands_Hphi = -99;
|
| 343 |
|
|
//FIXME
|
| 344 |
|
|
// float SimBs_dr = -99, SimBs_dEta= -99, SimBs_dPhi = -99, SimBs_Hmass = -99, SimBs_Hpt = -99, SimBs_Heta = -99, SimBs_Hphi = -99;
|
| 345 |
|
|
|
| 346 |
|
|
int Vtype,nSvs,nSimBs,numJets,numBJets,eventFlav;
|
| 347 |
|
|
// bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
|
| 348 |
|
|
bool triggerFlags[500],hbhe;
|
| 349 |
arizzi |
1.17 |
float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.;
|
| 350 |
arizzi |
1.1 |
// ----------------------------------------------------------------------
|
| 351 |
|
|
// First Part:
|
| 352 |
|
|
//
|
| 353 |
|
|
// * enable the AutoLibraryLoader
|
| 354 |
|
|
// * book the histograms of interest
|
| 355 |
|
|
// * open the input file
|
| 356 |
|
|
// ----------------------------------------------------------------------
|
| 357 |
|
|
|
| 358 |
|
|
// load framework libraries
|
| 359 |
bortigno |
1.30 |
gSystem->Load("libFWCoreFWLite");
|
| 360 |
arizzi |
1.1 |
gSystem->Load("libDataFormatsFWLite");
|
| 361 |
|
|
AutoLibraryLoader::enable();
|
| 362 |
|
|
|
| 363 |
|
|
// parse arguments
|
| 364 |
|
|
if ( argc < 2 ) {
|
| 365 |
|
|
return 0;
|
| 366 |
|
|
}
|
| 367 |
|
|
|
| 368 |
arizzi |
1.4 |
std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
|
| 369 |
|
|
std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
|
| 370 |
|
|
|
| 371 |
arizzi |
1.1 |
// get the python configuration
|
| 372 |
|
|
PythonProcessDesc builder(argv[1]);
|
| 373 |
|
|
const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
|
| 374 |
|
|
const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
|
| 375 |
|
|
const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
|
| 376 |
bortigno |
1.30 |
|
| 377 |
nmohr |
1.7 |
std::vector<edm::LuminosityBlockRange> jsonVector;
|
| 378 |
|
|
if ( in.exists("lumisToProcess") )
|
| 379 |
bortigno |
1.30 |
{
|
| 380 |
|
|
std::vector<edm::LuminosityBlockRange> const & lumisTemp =
|
| 381 |
|
|
in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
|
| 382 |
nmohr |
1.7 |
jsonVector.resize( lumisTemp.size() );
|
| 383 |
|
|
copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
|
| 384 |
bortigno |
1.30 |
}
|
| 385 |
arizzi |
1.1 |
|
| 386 |
|
|
// now get each parameter
|
| 387 |
|
|
int maxEvents_( in.getParameter<int>("maxEvents") );
|
| 388 |
nmohr |
1.29 |
int skipEvents_( in.getParameter<int>("skipEvents") );
|
| 389 |
arizzi |
1.1 |
unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
|
| 390 |
|
|
std::string outputFile_( out.getParameter<std::string>("fileName" ) );
|
| 391 |
|
|
std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
|
| 392 |
arizzi |
1.6 |
double btagThr = ana.getParameter<double>("bJetCountThreshold" );
|
| 393 |
arizzi |
1.4 |
bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
|
| 394 |
arizzi |
1.17 |
bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
|
| 395 |
|
|
bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
|
| 396 |
nmohr |
1.31 |
HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ);
|
| 397 |
|
|
HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW );
|
| 398 |
arizzi |
1.4 |
|
| 399 |
arizzi |
1.16 |
TriggerWeight triggerWeight(ana);
|
| 400 |
arizzi |
1.17 |
BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
|
| 401 |
|
|
BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
|
| 402 |
arizzi |
1.1 |
|
| 403 |
|
|
std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
|
| 404 |
bortigno |
1.30 |
// std::string inputFile( in.getParameter<std::string> ("fileName") );
|
| 405 |
arizzi |
1.1 |
|
| 406 |
arizzi |
1.3 |
std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
|
| 407 |
|
|
std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
|
| 408 |
arizzi |
1.25 |
|
| 409 |
arizzi |
1.1 |
bool isMC_( ana.getParameter<bool>("isMC") );
|
| 410 |
|
|
TriggerReader trigger(isMC_);
|
| 411 |
arizzi |
1.18 |
TriggerReader patFilters(false);
|
| 412 |
arizzi |
1.25 |
|
| 413 |
|
|
edm::LumiReWeighting lumiWeights;
|
| 414 |
|
|
if(isMC_)
|
| 415 |
bortigno |
1.30 |
{
|
| 416 |
|
|
// lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
|
| 417 |
|
|
}
|
| 418 |
arizzi |
1.1 |
|
| 419 |
bortigno |
1.30 |
// TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
|
| 420 |
|
|
// TH1F * pu = new TH1F("pileup","",-0.5,24.5,25);
|
| 421 |
|
|
TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
|
| 422 |
|
|
TH1F * count = new TH1F("Count","Count", 1,0,2 );
|
| 423 |
|
|
TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
|
| 424 |
arizzi |
1.1 |
_outTree = new TTree("tree", "myTree");
|
| 425 |
|
|
|
| 426 |
|
|
_outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F");
|
| 427 |
|
|
_outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
|
| 428 |
|
|
_outTree->Branch("nhJets" , &nhJets , "nhJets/I");
|
| 429 |
|
|
_outTree->Branch("naJets" , &naJets , "naJets/I");
|
| 430 |
arizzi |
1.2 |
|
| 431 |
|
|
_outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
|
| 432 |
|
|
_outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
|
| 433 |
|
|
_outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
|
| 434 |
|
|
_outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
|
| 435 |
|
|
_outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
|
| 436 |
|
|
_outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
|
| 437 |
|
|
_outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
|
| 438 |
|
|
_outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
|
| 439 |
|
|
_outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
|
| 440 |
|
|
_outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
|
| 441 |
|
|
_outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
|
| 442 |
|
|
_outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
|
| 443 |
|
|
_outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
|
| 444 |
|
|
_outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
|
| 445 |
|
|
_outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
|
| 446 |
|
|
_outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
|
| 447 |
|
|
_outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
|
| 448 |
arizzi |
1.24 |
_outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[naJets]/F");
|
| 449 |
|
|
_outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[naJets]/F");
|
| 450 |
|
|
_outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[naJets]/F");
|
| 451 |
arizzi |
1.2 |
|
| 452 |
|
|
_outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
|
| 453 |
|
|
_outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
|
| 454 |
|
|
_outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
|
| 455 |
|
|
_outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
|
| 456 |
|
|
_outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
|
| 457 |
|
|
_outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
|
| 458 |
|
|
_outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
|
| 459 |
|
|
_outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
|
| 460 |
|
|
_outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
|
| 461 |
|
|
_outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
|
| 462 |
|
|
_outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
|
| 463 |
|
|
_outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
|
| 464 |
|
|
_outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
|
| 465 |
|
|
_outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
|
| 466 |
|
|
_outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
|
| 467 |
|
|
_outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
|
| 468 |
|
|
_outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
|
| 469 |
arizzi |
1.24 |
_outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
|
| 470 |
|
|
_outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
|
| 471 |
|
|
_outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
|
| 472 |
arizzi |
1.1 |
|
| 473 |
|
|
_outTree->Branch("jjdr" , &jjdr , "jjdr/F" );
|
| 474 |
|
|
_outTree->Branch("jjdPhi" , &jjdPhi , "jjdPhi/F" );
|
| 475 |
arizzi |
1.6 |
_outTree->Branch("jjdEta" , &jjdEta , "jjdEta/F" );
|
| 476 |
arizzi |
1.1 |
_outTree->Branch("numJets" , &numJets , "numJets/I" );
|
| 477 |
|
|
_outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
|
| 478 |
|
|
_outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
|
| 479 |
|
|
_outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
|
| 480 |
|
|
_outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
|
| 481 |
|
|
_outTree->Branch("genZpt" , &genZpt , "genZpt/F");
|
| 482 |
|
|
_outTree->Branch("genWpt" , &genWpt , "genWpt/F");
|
| 483 |
|
|
_outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
|
| 484 |
arizzi |
1.22 |
_outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
|
| 485 |
arizzi |
1.24 |
_outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
|
| 486 |
|
|
_outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
|
| 487 |
|
|
_outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
|
| 488 |
|
|
|
| 489 |
|
|
|
| 490 |
arizzi |
1.1 |
_outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
|
| 491 |
|
|
_outTree->Branch("PUweight", &PUweight , "PUweight/F");
|
| 492 |
arizzi |
1.3 |
_outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
|
| 493 |
arizzi |
1.1 |
|
| 494 |
|
|
|
| 495 |
|
|
|
| 496 |
|
|
|
| 497 |
|
|
_outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
|
| 498 |
|
|
_outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
|
| 499 |
arizzi |
1.19 |
_outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
|
| 500 |
arizzi |
1.1 |
_outTree->Branch("VMt" , &VMt , "VMt/F" );
|
| 501 |
|
|
|
| 502 |
arizzi |
1.9 |
_outTree->Branch("nvlep" , &nvlep , "nvlep/I");
|
| 503 |
|
|
_outTree->Branch("nalep" , &nalep , "nalep/I");
|
| 504 |
arizzi |
1.3 |
|
| 505 |
arizzi |
1.9 |
_outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
|
| 506 |
|
|
_outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
|
| 507 |
|
|
_outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
|
| 508 |
|
|
_outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
|
| 509 |
|
|
_outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
|
| 510 |
|
|
_outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
|
| 511 |
|
|
_outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
|
| 512 |
|
|
_outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
|
| 513 |
|
|
_outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
|
| 514 |
|
|
_outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
|
| 515 |
|
|
_outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
|
| 516 |
|
|
_outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
|
| 517 |
|
|
_outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
|
| 518 |
arizzi |
1.27 |
_outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
|
| 519 |
|
|
_outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
|
| 520 |
arizzi |
1.9 |
|
| 521 |
|
|
_outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
|
| 522 |
|
|
_outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
|
| 523 |
|
|
_outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
|
| 524 |
|
|
_outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
|
| 525 |
|
|
_outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
|
| 526 |
|
|
_outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
|
| 527 |
|
|
_outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
|
| 528 |
|
|
_outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
|
| 529 |
|
|
_outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
|
| 530 |
|
|
_outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
|
| 531 |
|
|
_outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
|
| 532 |
|
|
_outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
|
| 533 |
|
|
_outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
|
| 534 |
arizzi |
1.27 |
_outTree->Branch("aLepton_id80",aLeptons.id80 ,"id[nalep]/F");
|
| 535 |
|
|
_outTree->Branch("aLepton_id95",aLeptons.id95 ,"id[nalep]/F");
|
| 536 |
arizzi |
1.9 |
|
| 537 |
arizzi |
1.8 |
_outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
|
| 538 |
arizzi |
1.1 |
|
| 539 |
bortigno |
1.30 |
//IVF
|
| 540 |
|
|
_outTree->Branch("nSvs",&nSvs ,"nSvs/I");
|
| 541 |
|
|
_outTree->Branch("BCands_massBCand", &IVF.massBcand,"massBcand[nSvs]/F");
|
| 542 |
|
|
_outTree->Branch("BCands_massSv", &IVF.massSv,"massSv[nSvs]/F");
|
| 543 |
|
|
_outTree->Branch("BCands_pt", &IVF.pt,"pt[nSvs]/F");
|
| 544 |
|
|
_outTree->Branch("BCands_eta", &IVF.eta,"eta[nSvs]/F");
|
| 545 |
|
|
_outTree->Branch("BCands_phi", &IVF.phi,"phi[nSvs]/F");
|
| 546 |
|
|
_outTree->Branch("BCands_dist3D", &IVF.dist3D,"dist3D[nSvs]/F");
|
| 547 |
|
|
_outTree->Branch("BCands_dist2D", &IVF.dist2D,"dist2D[nSvs]/F");
|
| 548 |
|
|
_outTree->Branch("BCands_distSim2D", &IVF.distSig2D,"distSig2D[nSvs]/F");
|
| 549 |
|
|
_outTree->Branch("BCands_distSig3D", &IVF.distSig3D,"distSig3D[nSvs]/F");
|
| 550 |
|
|
_outTree->Branch("BCands_dist3D_norm", &IVF.dist3D_norm,"dist3D_norm[nSvs]/F");
|
| 551 |
|
|
//IVF higgs candidate
|
| 552 |
|
|
_outTree->Branch("BCands_dr", &BCands_dr,"BCands_dr/F");
|
| 553 |
|
|
_outTree->Branch("BCands_dPhi", &BCands_dPhi,"BCands_dPhi/F");
|
| 554 |
|
|
_outTree->Branch("BCands_dEta", &BCands_dEta,"BCands_dEta/F");
|
| 555 |
|
|
_outTree->Branch("BCands_Hmass", &BCands_Hmass,"BCands_Hmass/F");
|
| 556 |
|
|
_outTree->Branch("BCands_Hpt", &BCands_Hpt,"BCands_Hpt/F");
|
| 557 |
|
|
_outTree->Branch("BCands_Heta", &BCands_Heta,"BCands_Heta/F");
|
| 558 |
|
|
_outTree->Branch("BCands_Hphi", &BCands_Hphi,"BCands_Hphi/F");
|
| 559 |
|
|
|
| 560 |
|
|
|
| 561 |
|
|
//FIXME : need to update the EDMNtuple with BHadron infos
|
| 562 |
|
|
// //SimBHadron
|
| 563 |
|
|
// _outTree->Branch("nSimBs",&nSimBs ,"nSimBs/I");
|
| 564 |
|
|
// _outTree->Branch("SimBs_mass", &SimBs.mass,"mass[nSvs]/F");
|
| 565 |
|
|
// _outTree->Branch("SimBs_pt", &SimBs.pt,"pt[nSvs]/F");
|
| 566 |
|
|
// _outTree->Branch("SimBs_eta", &SimBs.eta,"eta[nSvs]/F");
|
| 567 |
|
|
// _outTree->Branch("SimBs_phi", &SimBs.phi,"phi[nSvs]/F");
|
| 568 |
|
|
// _outTree->Branch("SimBs_vtx_x", &SimBs.vtx_x,"vtx_x[nSvs]/F");
|
| 569 |
|
|
// _outTree->Branch("SimBs_vtx_y", &SimBs.vtx_y,"vtx_y[nSvs]/F");
|
| 570 |
|
|
// _outTree->Branch("SimBs_vtx_z", &SimBs.vtx_z,"vtx_z[nSvs]/F");
|
| 571 |
|
|
// _outTree->Branch("SimBs_pdgId", &SimBs.pdgId,"pdgId[nSvs]/F");
|
| 572 |
|
|
// _outTree->Branch("SimBs_status", &SimBs.status,"status[nSvs]/F");
|
| 573 |
|
|
// //SimBHadron Higgs Candidate
|
| 574 |
|
|
// _outTree->Branch("SimBs_dr", &SimBs_dr,"SimBs_dr/F");
|
| 575 |
|
|
// _outTree->Branch("SimBs_dPhi", &SimBs_dPhi,"SimBs_dPhi/F");
|
| 576 |
|
|
// _outTree->Branch("SimBs_dEta", &SimBs_dEta,"SimBs_dEta/F");
|
| 577 |
|
|
// _outTree->Branch("SimBs_Hmass", &SimBs_Hmass,"SimBs_Hmass/F");
|
| 578 |
|
|
// _outTree->Branch("SimBs_Hpt", &SimBs_Hpt,"SimBs_Hpt/F");
|
| 579 |
|
|
// _outTree->Branch("SimBs_Heta", &SimBs_Heta,"SimBs_Heta/F");
|
| 580 |
|
|
// _outTree->Branch("SimBs_Hphi", &SimBs_Hphi,"SimBs_Hphi/F");
|
| 581 |
|
|
|
| 582 |
|
|
|
| 583 |
arizzi |
1.1 |
_outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
|
| 584 |
|
|
_outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
|
| 585 |
|
|
_outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
|
| 586 |
|
|
_outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
|
| 587 |
|
|
|
| 588 |
bortigno |
1.30 |
std::stringstream s;
|
| 589 |
|
|
s << "triggerFlags[" << triggers.size() << "]/b";
|
| 590 |
|
|
_outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
|
| 591 |
nmohr |
1.7 |
|
| 592 |
arizzi |
1.17 |
|
| 593 |
bortigno |
1.30 |
_outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
|
| 594 |
|
|
_outTree->Branch("hbhe" , &hbhe , "hbhe/b");
|
| 595 |
|
|
_outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
|
| 596 |
|
|
_outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
|
| 597 |
|
|
_outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
|
| 598 |
arizzi |
1.1 |
|
| 599 |
bortigno |
1.30 |
int ievt=0;
|
| 600 |
|
|
int totalcount=0;
|
| 601 |
arizzi |
1.1 |
|
| 602 |
bortigno |
1.30 |
// TFile* inFile = new TFile(inputFile.c_str(), "read");
|
| 603 |
arizzi |
1.1 |
for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
|
| 604 |
|
|
std::cout << iFile << std::endl;
|
| 605 |
|
|
TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
|
| 606 |
|
|
if(inFile==0) continue;
|
| 607 |
|
|
|
| 608 |
bortigno |
1.30 |
// loop the events
|
| 609 |
arizzi |
1.1 |
|
| 610 |
|
|
fwlite::Event ev(inFile);
|
| 611 |
nmohr |
1.29 |
for(ev.toBegin(); !ev.atEnd() ; ++ev, ++ievt)
|
| 612 |
arizzi |
1.1 |
{
|
| 613 |
nmohr |
1.29 |
if (ievt <= skipEvents_) continue;
|
| 614 |
|
|
if (maxEvents_ >= 0){
|
| 615 |
|
|
if (ievt > maxEvents_ + skipEvents_) break;
|
| 616 |
|
|
};
|
| 617 |
arizzi |
1.16 |
count->Fill(1.);
|
| 618 |
arizzi |
1.1 |
|
| 619 |
arizzi |
1.25 |
fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
|
| 620 |
|
|
vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 621 |
|
|
const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
|
| 622 |
|
|
PUweight=1.;
|
| 623 |
arizzi |
1.1 |
if(isMC_){
|
| 624 |
bortigno |
1.30 |
|
| 625 |
arizzi |
1.1 |
// PU weights
|
| 626 |
bortigno |
1.30 |
std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
|
| 627 |
|
|
PUweight = 1;//lumiWeights.weight( puit->second /3. );
|
| 628 |
|
|
}
|
| 629 |
|
|
countWithPU->Fill(PUweight);
|
| 630 |
nmohr |
1.7 |
|
| 631 |
bortigno |
1.30 |
//Write event info
|
| 632 |
|
|
EVENT.run = ev.id().run();
|
| 633 |
|
|
EVENT.lumi = ev.id().luminosityBlock();
|
| 634 |
|
|
EVENT.event = ev.id().event();
|
| 635 |
|
|
EVENT.json = jsonContainsEvent (jsonVector, ev);
|
| 636 |
|
|
|
| 637 |
|
|
// secondary vtxs
|
| 638 |
|
|
fwlite::Handle<std::vector<reco::Vertex> > SVC;
|
| 639 |
|
|
SVC.getByLabel(ev,"selectedVertices");
|
| 640 |
|
|
const std::vector<reco::Vertex> svc = *(SVC.product());
|
| 641 |
|
|
|
| 642 |
|
|
//FIXME : need to update EDM ntuple with BHadron infos
|
| 643 |
|
|
// // simBHadrons
|
| 644 |
|
|
// fwlite::Handle<SimBHadronCollection> SBHC;
|
| 645 |
|
|
// SBHC.getByLabel(ev, "bhadrons");
|
| 646 |
|
|
// const SimBHadronCollection sbhc = *(SBHC.product());
|
| 647 |
|
|
|
| 648 |
|
|
const std::vector<VHbbCandidate> * candZ ;
|
| 649 |
|
|
const std::vector<VHbbCandidate> * candW ;
|
| 650 |
|
|
const VHbbEvent * iEvent =0;
|
| 651 |
|
|
if(fromCandidate)
|
| 652 |
|
|
{
|
| 653 |
|
|
fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
|
| 654 |
|
|
vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
|
| 655 |
|
|
candZ = vhbbCandHandleZ.product();
|
| 656 |
|
|
|
| 657 |
|
|
fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
|
| 658 |
|
|
vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
|
| 659 |
|
|
candW = vhbbCandHandle.product();
|
| 660 |
|
|
}
|
| 661 |
|
|
else
|
| 662 |
|
|
{
|
| 663 |
|
|
candZlocal->clear();
|
| 664 |
|
|
candWlocal->clear();
|
| 665 |
|
|
fwlite::Handle< VHbbEvent > vhbbHandle;
|
| 666 |
|
|
vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 667 |
|
|
iEvent = vhbbHandle.product();
|
| 668 |
|
|
algoZ->run(vhbbHandle.product(),*candZlocal);
|
| 669 |
|
|
algoW->run(vhbbHandle.product(),*candWlocal);
|
| 670 |
|
|
candZ= candZlocal;
|
| 671 |
|
|
candW= candWlocal;
|
| 672 |
|
|
|
| 673 |
|
|
/* for(size_t m=0;m<iEvent->muInfo.size();m++)
|
| 674 |
|
|
{
|
| 675 |
|
|
|
| 676 |
|
|
if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
|
| 677 |
|
|
fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
|
| 678 |
|
|
{
|
| 679 |
|
|
std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
|
| 680 |
|
|
}
|
| 681 |
|
|
}
|
| 682 |
|
|
*/
|
| 683 |
arizzi |
1.4 |
|
| 684 |
bortigno |
1.30 |
}
|
| 685 |
arizzi |
1.1 |
|
| 686 |
arizzi |
1.3 |
const std::vector<VHbbCandidate> * cand = candZ;
|
| 687 |
arizzi |
1.1 |
|
| 688 |
|
|
|
| 689 |
|
|
|
| 690 |
bortigno |
1.30 |
/* fwlite::Handle< VHbbEvent > vhbbHandle;
|
| 691 |
|
|
vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 692 |
|
|
const VHbbEvent iEvent = *vhbbHandle.product();
|
| 693 |
|
|
*/
|
| 694 |
arizzi |
1.1 |
|
| 695 |
bortigno |
1.30 |
// std::clog << "Filling tree "<< std::endl;
|
| 696 |
|
|
bool isW=false;
|
| 697 |
arizzi |
1.21 |
|
| 698 |
arizzi |
1.17 |
|
| 699 |
bortigno |
1.30 |
if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
|
| 700 |
|
|
if(cand->size() > 1 )
|
| 701 |
arizzi |
1.3 |
{
|
| 702 |
bortigno |
1.30 |
std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
|
| 703 |
arizzi |
1.3 |
}
|
| 704 |
bortigno |
1.30 |
if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
|
| 705 |
|
|
if(cand->size() == 0)
|
| 706 |
arizzi |
1.3 |
{
|
| 707 |
bortigno |
1.30 |
// std::cout << "W event loss due to tigther cuts" << std::endl;
|
| 708 |
arizzi |
1.3 |
continue;
|
| 709 |
|
|
}
|
| 710 |
bortigno |
1.30 |
const VHbbCandidate & vhCand = cand->at(0);
|
| 711 |
|
|
patFilters.setEvent(&ev,"VH");
|
| 712 |
|
|
hbhe = patFilters.accept("hbhe");
|
| 713 |
arizzi |
1.18 |
|
| 714 |
bortigno |
1.30 |
trigger.setEvent(&ev);
|
| 715 |
|
|
for(size_t j=0;j < triggers.size();j++)
|
| 716 |
arizzi |
1.4 |
triggerFlags[j]=trigger.accept(triggers[j]);
|
| 717 |
arizzi |
1.3 |
|
| 718 |
bortigno |
1.30 |
eventFlav=0;
|
| 719 |
|
|
if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
|
| 720 |
|
|
else if(aux.mcC.size() > 0) eventFlav=4;
|
| 721 |
arizzi |
1.8 |
|
| 722 |
arizzi |
1.3 |
|
| 723 |
bortigno |
1.30 |
H.mass = vhCand.H.p4.M();
|
| 724 |
|
|
H.pt = vhCand.H.p4.Pt();
|
| 725 |
|
|
H.eta = vhCand.H.p4.Eta();
|
| 726 |
|
|
H.phi = vhCand.H.p4.Phi();
|
| 727 |
|
|
V.mass = vhCand.V.p4.M();
|
| 728 |
|
|
V.pt = vhCand.V.p4.Pt();
|
| 729 |
|
|
V.eta = vhCand.V.p4.Eta();
|
| 730 |
|
|
V.phi = vhCand.V.p4.Phi();
|
| 731 |
|
|
nhJets=2;
|
| 732 |
|
|
hJets.set(vhCand.H.jets[0],0);
|
| 733 |
|
|
hJets.set(vhCand.H.jets[1],1);
|
| 734 |
|
|
aJets.reset();
|
| 735 |
|
|
naJets=vhCand.additionalJets.size();
|
| 736 |
|
|
numBJets=0;
|
| 737 |
|
|
if(vhCand.H.jets[0].csv> btagThr) numBJets++;
|
| 738 |
|
|
if(vhCand.H.jets[1].csv> btagThr) numBJets++;
|
| 739 |
|
|
for( int j=0; j < naJets && j < MAXJ; j++ )
|
| 740 |
arizzi |
1.6 |
{
|
| 741 |
bortigno |
1.30 |
aJets.set(vhCand.additionalJets[j],j);
|
| 742 |
|
|
if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
|
| 743 |
arizzi |
1.6 |
}
|
| 744 |
bortigno |
1.30 |
numJets = vhCand.additionalJets.size()+2;
|
| 745 |
|
|
jjdr = deltaR(vhCand.H.jets[0].p4.Eta(),vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Eta(),vhCand.H.jets[1].p4.Phi());
|
| 746 |
|
|
jjdPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
|
| 747 |
|
|
jjdEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
|
| 748 |
|
|
HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
|
| 749 |
|
|
HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
|
| 750 |
|
|
VMt = vhCand.Mt() ;
|
| 751 |
|
|
deltaPullAngle = vhCand.H.deltaTheta;
|
| 752 |
|
|
|
| 753 |
|
|
hJets.cosTheta[0]= vhCand.H.helicities[0];
|
| 754 |
|
|
hJets.cosTheta[1]= vhCand.H.helicities[1];
|
| 755 |
arizzi |
1.24 |
|
| 756 |
bortigno |
1.30 |
MET.et = vhCand.V.mets.at(0).p4.Pt();
|
| 757 |
|
|
MET.phi = vhCand.V.mets.at(0).p4.Phi();
|
| 758 |
|
|
MET.sumet = vhCand.V.mets.at(0).sumEt;
|
| 759 |
|
|
MET.sig = vhCand.V.mets.at(0).metSig;
|
| 760 |
arizzi |
1.24 |
|
| 761 |
|
|
if(!fromCandidate) {
|
| 762 |
|
|
MHT.mht = iEvent->mht.p4.Pt();
|
| 763 |
|
|
MHT.phi = iEvent->mht.p4.Phi();
|
| 764 |
|
|
MHT.ht = iEvent->mht.sumEt;
|
| 765 |
|
|
MHT.sig = iEvent->mht.metSig;
|
| 766 |
bortigno |
1.30 |
}
|
| 767 |
|
|
|
| 768 |
|
|
Vtype = vhCand.candidateType;
|
| 769 |
arizzi |
1.24 |
|
| 770 |
|
|
|
| 771 |
bortigno |
1.30 |
//Secondary Vertices
|
| 772 |
|
|
IVF.reset();
|
| 773 |
|
|
nSvs = svc.size();
|
| 774 |
|
|
const TVector3 recoPv = aux.pvInfo.firstPVInPT2;
|
| 775 |
|
|
const math::XYZPoint myPv(recoPv);
|
| 776 |
|
|
//look here for Matrix filling info http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/SMatrixDoc.html
|
| 777 |
|
|
std::vector<double> fillMatrix(6);
|
| 778 |
|
|
for (int i = 0; i<6; ++i) fillMatrix[i] = 0.;
|
| 779 |
|
|
fillMatrix[0] = TMath::Power(0.002,2);
|
| 780 |
|
|
fillMatrix[2] = TMath::Power(0.002,2);
|
| 781 |
|
|
fillMatrix[5] = TMath::Power(0.002,2);
|
| 782 |
|
|
const ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > myFakeMatrixError(fillMatrix.begin(),fillMatrix.end());
|
| 783 |
|
|
const reco::Vertex recoVtxPv(myPv, myFakeMatrixError);
|
| 784 |
|
|
for( int j=0; j < nSvs && j < MAXB; ++j ) {
|
| 785 |
|
|
const GlobalVector flightDir = flightDirection(recoPv,svc[j]);
|
| 786 |
|
|
reco::SecondaryVertex recoSv(recoVtxPv, svc[j], flightDir ,true);
|
| 787 |
|
|
IVF.set( recoSv, recoPv ,j);
|
| 788 |
|
|
}
|
| 789 |
|
|
if(nSvs > 2){
|
| 790 |
|
|
TLorentzVector BCands_H1, BCands_H2, BCands_H;
|
| 791 |
|
|
BCands_H1.SetPtEtaPhiM(IVF.pt[0], IVF.eta[0], IVF.phi[0], IVF.massBcand[0]);
|
| 792 |
|
|
BCands_H2.SetPtEtaPhiM(IVF.pt[1], IVF.eta[1], IVF.phi[1], IVF.massBcand[1]);
|
| 793 |
|
|
BCands_H = BCands_H1 + BCands_H2;
|
| 794 |
|
|
BCands_dr = deltaR(IVF.eta[0], IVF.phi[0], IVF.eta[1], IVF.phi[1] );
|
| 795 |
|
|
BCands_dPhi = deltaPhi(IVF.phi[0], IVF.phi[1] );
|
| 796 |
|
|
BCands_dEta = TMath::Abs(IVF.eta[0] - IVF.eta[1] );
|
| 797 |
|
|
BCands_Hmass = BCands_H.M();
|
| 798 |
|
|
BCands_Hpt = BCands_H.Pt();
|
| 799 |
|
|
BCands_Heta = BCands_H.Eta();
|
| 800 |
|
|
BCands_Hphi = BCands_H.Phi();
|
| 801 |
|
|
}
|
| 802 |
|
|
|
| 803 |
|
|
//FIXME : need to update EDM ntuple with simBhadron info
|
| 804 |
|
|
// //SimBhadron
|
| 805 |
|
|
// SimBs.reset();
|
| 806 |
|
|
// nSimBs = sbhc.size();
|
| 807 |
|
|
// for( int j=0; j < nSimBs && j < MAXB; ++j )
|
| 808 |
|
|
// SimBs.set( sbhc.at(j), j);
|
| 809 |
|
|
// if(nSimBs > 2){
|
| 810 |
|
|
// TLorentzVector SimBs_H1, SimBs_H2, SimBs_H;
|
| 811 |
|
|
// SimBs_H1.SetPtEtaPhiM(SimBs.pt[0], SimBs.eta[0], SimBs.phi[0], SimBs.mass[0]);
|
| 812 |
|
|
// SimBs_H2.SetPtEtaPhiM(SimBs.pt[1], SimBs.eta[1], SimBs.phi[1], SimBs.mass[1]);
|
| 813 |
|
|
// SimBs_H = SimBs_H1 + SimBs_H2;
|
| 814 |
|
|
// SimBs_dr = deltaR(SimBs.eta[0], SimBs.phi[0], SimBs.eta[1], SimBs.phi[1] );
|
| 815 |
|
|
// SimBs_dPhi = deltaPhi(SimBs.phi[0], SimBs.phi[1] );
|
| 816 |
|
|
// SimBs_dEta = TMath::Abs(SimBs.eta[0] - SimBs.eta[1] );
|
| 817 |
|
|
// SimBs_Hmass = SimBs_H.M();
|
| 818 |
|
|
// SimBs_Hpt = SimBs_H.Pt();
|
| 819 |
|
|
// SimBs_Heta = SimBs_H.Eta();
|
| 820 |
|
|
// SimBs_Hphi = SimBs_H.Phi();
|
| 821 |
|
|
// }
|
| 822 |
|
|
|
| 823 |
|
|
|
| 824 |
|
|
//Loop on jets
|
| 825 |
|
|
double maxBtag=-99999;
|
| 826 |
|
|
minDeltaPhijetMET = 999;
|
| 827 |
|
|
TLorentzVector bJet;
|
| 828 |
|
|
std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
|
| 829 |
|
|
std::vector<float> jet30eta;
|
| 830 |
|
|
std::vector<float> jet30pt;
|
| 831 |
|
|
if(fromCandidate)
|
| 832 |
|
|
{
|
| 833 |
|
|
//Loop on Higgs Jets
|
| 834 |
|
|
for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
|
| 835 |
|
|
if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
|
| 836 |
|
|
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()) < minDeltaPhijetMET)
|
| 837 |
|
|
{
|
| 838 |
arizzi |
1.24 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi());
|
| 839 |
|
|
jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
|
| 840 |
bortigno |
1.30 |
}
|
| 841 |
|
|
btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
|
| 842 |
|
|
}
|
| 843 |
|
|
//Loop on Additional Jets
|
| 844 |
|
|
for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
|
| 845 |
|
|
if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
|
| 846 |
|
|
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()) < minDeltaPhijetMET)
|
| 847 |
|
|
{
|
| 848 |
arizzi |
1.24 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi());
|
| 849 |
|
|
jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
|
| 850 |
bortigno |
1.30 |
}
|
| 851 |
|
|
if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
|
| 852 |
|
|
{
|
| 853 |
|
|
if(vhCand.additionalJets[j].p4.Pt() > 20)
|
| 854 |
|
|
btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
|
| 855 |
|
|
}
|
| 856 |
|
|
}
|
| 857 |
|
|
}
|
| 858 |
|
|
else
|
| 859 |
|
|
{
|
| 860 |
|
|
for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
|
| 861 |
|
|
if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
|
| 862 |
|
|
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()) < minDeltaPhijetMET)
|
| 863 |
|
|
{
|
| 864 |
arizzi |
1.24 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi());
|
| 865 |
|
|
jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
|
| 866 |
bortigno |
1.30 |
}
|
| 867 |
|
|
if(iEvent->simpleJets2[j].p4.Pt() > 30)
|
| 868 |
|
|
{
|
| 869 |
|
|
jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
|
| 870 |
|
|
jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
|
| 871 |
|
|
}
|
| 872 |
arizzi |
1.24 |
|
| 873 |
bortigno |
1.30 |
//For events made with highest CSV, all jets in the event should be taken into account for "tagging" SF (anti tagging is a mess)
|
| 874 |
|
|
// because for example a light jet not used for the Higgs can have in reality a higher CSV due to SF > 1 and become a higgs jet
|
| 875 |
|
|
if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
|
| 876 |
|
|
{
|
| 877 |
|
|
if(iEvent->simpleJets2[j].p4.Pt() > 20)
|
| 878 |
|
|
btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
|
| 879 |
|
|
}
|
| 880 |
|
|
}
|
| 881 |
|
|
|
| 882 |
|
|
//if we use the highest pt pair, only the two higgs jet should be used to compute the SF because the other jets are excluded
|
| 883 |
|
|
// by a criteria (pt of the dijet) that is not btag SF dependent
|
| 884 |
|
|
if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
|
| 885 |
|
|
for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
|
| 886 |
|
|
}
|
| 887 |
arizzi |
1.24 |
|
| 888 |
bortigno |
1.30 |
}
|
| 889 |
|
|
vLeptons.reset();
|
| 890 |
|
|
weightTrig = 1.; // better to default to 1
|
| 891 |
|
|
TLorentzVector leptonForTop;
|
| 892 |
|
|
size_t firstAddMu=0;
|
| 893 |
|
|
size_t firstAddEle=0;
|
| 894 |
|
|
if(Vtype == VHbbCandidate::Zmumu ){
|
| 895 |
|
|
vLeptons.set(vhCand.V.muons[0],0,13);
|
| 896 |
|
|
vLeptons.set(vhCand.V.muons[1],1,13);
|
| 897 |
|
|
float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
|
| 898 |
|
|
float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
|
| 899 |
|
|
float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
|
| 900 |
|
|
float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
|
| 901 |
|
|
weightTrig = cweightID * cweightTrig;
|
| 902 |
|
|
nvlep=2;
|
| 903 |
|
|
firstAddMu=2;
|
| 904 |
|
|
}
|
| 905 |
|
|
if( Vtype == VHbbCandidate::Zee ){
|
| 906 |
|
|
vLeptons.set(vhCand.V.electrons[0],0,11);
|
| 907 |
|
|
vLeptons.set(vhCand.V.electrons[1],1,11);
|
| 908 |
|
|
nvlep=2;
|
| 909 |
|
|
firstAddEle=2;
|
| 910 |
|
|
std::vector<float> pt,eta;
|
| 911 |
|
|
pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
|
| 912 |
|
|
pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
|
| 913 |
|
|
weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
|
| 914 |
|
|
triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
|
| 915 |
|
|
weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
|
| 916 |
|
|
weightTrig = weightEleTrigElePart * weightEleRecoAndId;
|
| 917 |
arizzi |
1.24 |
|
| 918 |
|
|
|
| 919 |
|
|
|
| 920 |
bortigno |
1.30 |
}
|
| 921 |
|
|
if(Vtype == VHbbCandidate::Wmun ){
|
| 922 |
|
|
leptonForTop=vhCand.V.muons[0].p4;
|
| 923 |
|
|
vLeptons.set(vhCand.V.muons[0],0,13);
|
| 924 |
|
|
float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
|
| 925 |
|
|
float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
|
| 926 |
|
|
float cweightTrig = weightTrig1;
|
| 927 |
|
|
weightTrig = cweightID * cweightTrig;
|
| 928 |
|
|
nvlep=1;
|
| 929 |
|
|
firstAddMu=1;
|
| 930 |
|
|
}
|
| 931 |
|
|
if( Vtype == VHbbCandidate::Wen ){
|
| 932 |
|
|
leptonForTop=vhCand.V.electrons[0].p4;
|
| 933 |
|
|
vLeptons.set(vhCand.V.electrons[0],0,11);
|
| 934 |
|
|
nvlep=1;
|
| 935 |
|
|
firstAddEle=1;
|
| 936 |
|
|
float weightMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
|
| 937 |
|
|
float weightV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
|
| 938 |
|
|
weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
|
| 939 |
|
|
weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30eta,jet30pt)*triggerWeight.scalePFMHTEle(MET.et);
|
| 940 |
|
|
weightEleTrigElePart= weightV4; //this is for debugging only, checking only the V4 part
|
| 941 |
|
|
|
| 942 |
|
|
weightMay*=weightEleRecoAndId;
|
| 943 |
|
|
weightV4*=weightEleRecoAndId;
|
| 944 |
|
|
weightV4*=weightEleTrigJetMETPart;
|
| 945 |
|
|
weightTrig = weightMay * 0.187 + weightV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
|
| 946 |
|
|
}
|
| 947 |
|
|
if( Vtype == VHbbCandidate::Znn ){
|
| 948 |
|
|
nvlep=0;
|
| 949 |
|
|
float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
|
| 950 |
|
|
weightTrig = weightTrig1;
|
| 951 |
|
|
weightTrigMET = weightTrig1;
|
| 952 |
|
|
}
|
| 953 |
arizzi |
1.9 |
|
| 954 |
bortigno |
1.30 |
aLeptons.reset();
|
| 955 |
|
|
nalep=0;
|
| 956 |
|
|
if(fromCandidate)
|
| 957 |
arizzi |
1.21 |
{
|
| 958 |
|
|
for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
|
| 959 |
|
|
for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
|
| 960 |
|
|
}
|
| 961 |
bortigno |
1.30 |
else
|
| 962 |
arizzi |
1.21 |
{
|
| 963 |
bortigno |
1.30 |
for(size_t j=0;j< iEvent->muInfo.size();j++)
|
| 964 |
|
|
{
|
| 965 |
arizzi |
1.21 |
if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
|
| 966 |
bortigno |
1.30 |
aLeptons.set(iEvent->muInfo[j],nalep++,13);
|
| 967 |
|
|
}
|
| 968 |
|
|
for(size_t j=0;j< iEvent->eleInfo.size();j++)
|
| 969 |
|
|
{
|
| 970 |
arizzi |
1.21 |
if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
|
| 971 |
bortigno |
1.30 |
aLeptons.set(iEvent->eleInfo[j],nalep++,11);
|
| 972 |
|
|
}
|
| 973 |
arizzi |
1.8 |
|
| 974 |
arizzi |
1.21 |
}
|
| 975 |
arizzi |
1.8 |
|
| 976 |
arizzi |
1.24 |
|
| 977 |
bortigno |
1.30 |
if(isMC_)
|
| 978 |
|
|
{
|
| 979 |
|
|
//std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
|
| 980 |
|
|
if ( btagJetInfos.size()< 10)
|
| 981 |
|
|
{
|
| 982 |
|
|
btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
|
| 983 |
|
|
btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
|
| 984 |
|
|
btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
|
| 985 |
|
|
}
|
| 986 |
|
|
else
|
| 987 |
|
|
{
|
| 988 |
|
|
std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
|
| 989 |
|
|
//TODO: revert to random throw for this cases
|
| 990 |
|
|
btag1T2CSF = 1.;
|
| 991 |
|
|
btag2TSF = 1.;
|
| 992 |
|
|
btag1TSF = 1.;
|
| 993 |
|
|
}
|
| 994 |
|
|
}
|
| 995 |
arizzi |
1.17 |
|
| 996 |
bortigno |
1.30 |
if(maxBtag > -99999)
|
| 997 |
arizzi |
1.8 |
{
|
| 998 |
bortigno |
1.30 |
TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
|
| 999 |
|
|
top.mass = topQuark.p4.M();
|
| 1000 |
|
|
top.pt = topQuark.p4.Pt();
|
| 1001 |
|
|
top.wMass = topQuark.p4W.M();
|
| 1002 |
arizzi |
1.8 |
} else {
|
| 1003 |
bortigno |
1.30 |
top.mass = -99;
|
| 1004 |
|
|
top.pt = -99;
|
| 1005 |
|
|
top.wMass = -99;
|
| 1006 |
|
|
}
|
| 1007 |
arizzi |
1.8 |
|
| 1008 |
|
|
|
| 1009 |
|
|
|
| 1010 |
bortigno |
1.30 |
//FIXME: too much warnings... figure out why
|
| 1011 |
|
|
// gendrcc=aux.genCCDeltaR();
|
| 1012 |
|
|
// gendrbb=aux.genBBDeltaR();
|
| 1013 |
|
|
genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
|
| 1014 |
|
|
genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
|
| 1015 |
arizzi |
1.24 |
|
| 1016 |
arizzi |
1.25 |
|
| 1017 |
arizzi |
1.32 |
/// Compute pull angle from AK7
|
| 1018 |
arizzi |
1.25 |
if(!fromCandidate){
|
| 1019 |
arizzi |
1.32 |
std::vector<VHbbEvent::SimpleJet> ak7wrt1(iEvent->simpleJets3);
|
| 1020 |
|
|
std::vector<VHbbEvent::SimpleJet> ak7wrt2(iEvent->simpleJets3);
|
| 1021 |
|
|
if(ak7wrt1.size() > 1){
|
| 1022 |
arizzi |
1.24 |
CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
|
| 1023 |
|
|
CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
|
| 1024 |
arizzi |
1.32 |
std::sort( ak7wrt1.begin(),ak7wrt1.end(),deltaRComparatorJ1 );
|
| 1025 |
|
|
std::sort( ak7wrt2.begin(),ak7wrt2.end(),deltaRComparatorJ2 );
|
| 1026 |
arizzi |
1.24 |
std::vector<VHbbEvent::SimpleJet> ak7_matched;
|
| 1027 |
arizzi |
1.32 |
// if the matched are different save them
|
| 1028 |
|
|
if(ak7wrt1[0].p4.DeltaR(ak7wrt2[0].p4) > 0.1) {
|
| 1029 |
|
|
ak7_matched.push_back(ak7wrt1[0]);
|
| 1030 |
|
|
ak7_matched.push_back(ak7wrt2[0]);
|
| 1031 |
|
|
}
|
| 1032 |
|
|
// else look at the second best
|
| 1033 |
|
|
else{
|
| 1034 |
|
|
// ak7wrt1 is best
|
| 1035 |
|
|
if( ak7wrt1[1].p4.DeltaR(vhCand.H.jets[0].p4) < ak7wrt2[1].p4.DeltaR(vhCand.H.jets[1].p4))
|
| 1036 |
|
|
{
|
| 1037 |
|
|
ak7_matched.push_back(ak7wrt1[1]);
|
| 1038 |
|
|
ak7_matched.push_back(ak7wrt2[0]);
|
| 1039 |
|
|
}
|
| 1040 |
|
|
else
|
| 1041 |
|
|
{
|
| 1042 |
|
|
ak7_matched.push_back(ak7wrt1[0]);
|
| 1043 |
|
|
ak7_matched.push_back(ak7wrt2[1]);
|
| 1044 |
|
|
}
|
| 1045 |
|
|
}
|
| 1046 |
arizzi |
1.24 |
CompareJetPt ptComparator;
|
| 1047 |
|
|
std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
|
| 1048 |
|
|
if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
|
| 1049 |
arizzi |
1.32 |
and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5)
|
| 1050 |
|
|
{
|
| 1051 |
|
|
deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
|
| 1052 |
|
|
}
|
| 1053 |
arizzi |
1.24 |
}
|
| 1054 |
|
|
}
|
| 1055 |
arizzi |
1.3 |
|
| 1056 |
bortigno |
1.30 |
_outTree->Fill();
|
| 1057 |
arizzi |
1.4 |
|
| 1058 |
bortigno |
1.30 |
}// closed event loop
|
| 1059 |
arizzi |
1.1 |
|
| 1060 |
bortigno |
1.30 |
std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
|
| 1061 |
|
|
inFile->Close();
|
| 1062 |
|
|
// close input file
|
| 1063 |
|
|
} // loop on files
|
| 1064 |
arizzi |
1.1 |
|
| 1065 |
|
|
|
| 1066 |
bortigno |
1.30 |
std::cout << "Events: " << ievt <<std::endl;
|
| 1067 |
|
|
std::cout << "TotalCount: " << totalcount <<std::endl;
|
| 1068 |
arizzi |
1.1 |
|
| 1069 |
|
|
|
| 1070 |
|
|
|
| 1071 |
bortigno |
1.30 |
_outFile->cd();
|
| 1072 |
arizzi |
1.1 |
|
| 1073 |
bortigno |
1.30 |
_outTree->Write();
|
| 1074 |
|
|
_outFile->Write();
|
| 1075 |
|
|
_outFile->Close();
|
| 1076 |
|
|
return 0;
|
| 1077 |
arizzi |
1.1 |
}
|
| 1078 |
|
|
|
| 1079 |
|
|
|
