| 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 |
#include "VHbbAnalysis/VHbbDataFormats/interface/HbbCandidateFinderAlgo.h"
|
| 24 |
#include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc"
|
| 25 |
|
| 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 |
#include "VHbbAnalysis/VHbbDataFormats/interface/TopMassReco.h"
|
| 31 |
|
| 32 |
|
| 33 |
//Move class definition to Ntupler.h ?
|
| 34 |
//#include "VHbbAnalysis/VHbbDataFormats/interface/Ntupler.h"
|
| 35 |
|
| 36 |
#include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
|
| 37 |
#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
|
| 38 |
|
| 39 |
#include <sstream>
|
| 40 |
#include <string>
|
| 41 |
|
| 42 |
#define MAXJ 30
|
| 43 |
#define MAXL 10
|
| 44 |
|
| 45 |
struct CompareDeltaR {
|
| 46 |
CompareDeltaR(TLorentzVector p4dir_): p4dir(p4dir_) {}
|
| 47 |
bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
|
| 48 |
return j1.p4.DeltaR(p4dir) > j2.p4.DeltaR(p4dir);
|
| 49 |
}
|
| 50 |
TLorentzVector p4dir;
|
| 51 |
};
|
| 52 |
|
| 53 |
bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
|
| 54 |
const edm::EventBase &event)
|
| 55 |
{
|
| 56 |
// if the jsonVec is empty, then no JSON file was provided so all
|
| 57 |
// events should pass
|
| 58 |
if (jsonVec.empty())
|
| 59 |
{
|
| 60 |
return true;
|
| 61 |
}
|
| 62 |
bool (* funcPtr) (edm::LuminosityBlockRange const &,
|
| 63 |
edm::LuminosityBlockID const &) = &edm::contains;
|
| 64 |
edm::LuminosityBlockID lumiID (event.id().run(),
|
| 65 |
event.id().luminosityBlock());
|
| 66 |
std::vector< edm::LuminosityBlockRange >::const_iterator iter =
|
| 67 |
std::find_if (jsonVec.begin(), jsonVec.end(),
|
| 68 |
boost::bind(funcPtr, _1, lumiID) );
|
| 69 |
return jsonVec.end() != iter;
|
| 70 |
|
| 71 |
}
|
| 72 |
|
| 73 |
|
| 74 |
|
| 75 |
typedef struct
|
| 76 |
{
|
| 77 |
float mass; //MT in case of W
|
| 78 |
float pt;
|
| 79 |
float eta;
|
| 80 |
float phi;
|
| 81 |
} TrackInfo;
|
| 82 |
|
| 83 |
|
| 84 |
struct LeptonInfo
|
| 85 |
{
|
| 86 |
void reset()
|
| 87 |
{
|
| 88 |
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; }
|
| 89 |
}
|
| 90 |
|
| 91 |
template <class Input> void set(const Input & i, int j,int t)
|
| 92 |
{
|
| 93 |
type[j]=t;
|
| 94 |
pt[j]=i.p4.Pt();
|
| 95 |
mass[j]=i.p4.M();
|
| 96 |
eta[j]=i.p4.Eta();
|
| 97 |
phi[j]=i.p4.Phi();
|
| 98 |
aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
|
| 99 |
pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
|
| 100 |
photonIso[j]=i.pfPhoIso;
|
| 101 |
neutralHadIso[j]=i.pfNeuIso;
|
| 102 |
chargedHadIso[j]=i.pfChaIso;
|
| 103 |
setSpecific(i,j);
|
| 104 |
}
|
| 105 |
template <class Input> void setSpecific(const Input & i, int j)
|
| 106 |
{
|
| 107 |
}
|
| 108 |
|
| 109 |
|
| 110 |
|
| 111 |
|
| 112 |
float mass[MAXL]; //MT in case of W
|
| 113 |
float pt[MAXL];
|
| 114 |
float eta[MAXL];
|
| 115 |
float phi[MAXL];
|
| 116 |
float aodCombRelIso[MAXL];
|
| 117 |
float pfCombRelIso[MAXL];
|
| 118 |
float photonIso[MAXL];
|
| 119 |
float neutralHadIso[MAXL];
|
| 120 |
float chargedHadIso[MAXL];
|
| 121 |
float particleIso[MAXL];
|
| 122 |
float dxy[MAXL];
|
| 123 |
float dz[MAXL];
|
| 124 |
int type[MAXL];
|
| 125 |
float id80[MAXL];
|
| 126 |
float id95[MAXL];
|
| 127 |
};
|
| 128 |
|
| 129 |
template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
|
| 130 |
id80[j]=i.id80r;
|
| 131 |
id95[j]=i.id95r;
|
| 132 |
}
|
| 133 |
template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
|
| 134 |
dxy[j]=i.ipDb;
|
| 135 |
dz[j]=i.zPVPt;
|
| 136 |
}
|
| 137 |
|
| 138 |
typedef struct
|
| 139 |
{
|
| 140 |
float et;
|
| 141 |
float sumet;
|
| 142 |
float sig;
|
| 143 |
float phi;
|
| 144 |
} METInfo;
|
| 145 |
|
| 146 |
typedef struct
|
| 147 |
{
|
| 148 |
float mht;
|
| 149 |
float ht;
|
| 150 |
float sig;
|
| 151 |
float phi;
|
| 152 |
} MHTInfo;
|
| 153 |
|
| 154 |
typedef struct
|
| 155 |
{
|
| 156 |
float mass;
|
| 157 |
float pt;
|
| 158 |
float wMass;
|
| 159 |
} TopInfo;
|
| 160 |
|
| 161 |
typedef struct
|
| 162 |
{
|
| 163 |
int run;
|
| 164 |
int lumi;
|
| 165 |
int event;
|
| 166 |
int json;
|
| 167 |
} EventInfo;
|
| 168 |
|
| 169 |
typedef struct
|
| 170 |
{
|
| 171 |
void set(const VHbbEvent::SimpleJet & j, int i)
|
| 172 |
{
|
| 173 |
pt[i]=j.p4.Pt();
|
| 174 |
eta[i]=j.p4.Eta();
|
| 175 |
phi[i]=j.p4.Phi();
|
| 176 |
csv[i]=j.csv;
|
| 177 |
numTracksSV[i] = j.vtxMass;
|
| 178 |
vtxMass[i]= j.vtxMass;
|
| 179 |
vtx3dL[i] = j.vtx3dL;
|
| 180 |
vtx3deL[i] = j.vtx3deL;
|
| 181 |
chf[i]=j.chargedHadronEFraction;
|
| 182 |
nhf[i] =j.neutralHadronEFraction;
|
| 183 |
cef[i] =j.chargedEmEFraction;
|
| 184 |
nef[i] =j.neutralEmEFraction;
|
| 185 |
nconstituents[i] =j.nConstituents;
|
| 186 |
nch[i]=j.ntracks;
|
| 187 |
|
| 188 |
flavour[i]=j.flavour;
|
| 189 |
if(j.bestMCp4.Pt() > 0)
|
| 190 |
{
|
| 191 |
genPt[i]=j.bestMCp4.Pt();
|
| 192 |
genEta[i]=j.bestMCp4.Eta();
|
| 193 |
genPhi[i]=j.bestMCp4.Phi();
|
| 194 |
}
|
| 195 |
JECUnc[i]=j.jecunc;
|
| 196 |
|
| 197 |
}
|
| 198 |
void reset()
|
| 199 |
{
|
| 200 |
for(int i=0;i<MAXJ;i++) {
|
| 201 |
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;
|
| 202 |
}
|
| 203 |
}
|
| 204 |
float pt[MAXJ];
|
| 205 |
float eta[MAXJ];
|
| 206 |
float phi[MAXJ];
|
| 207 |
float csv[MAXJ];
|
| 208 |
float cosTheta[MAXJ];
|
| 209 |
int numTracksSV[MAXJ];
|
| 210 |
float chf[MAXJ];
|
| 211 |
float nhf[MAXJ];
|
| 212 |
float cef[MAXJ];
|
| 213 |
float nef[MAXJ];
|
| 214 |
float nch[MAXJ];
|
| 215 |
float nconstituents[MAXJ];
|
| 216 |
float flavour[MAXJ];
|
| 217 |
float genPt[MAXJ];
|
| 218 |
float genEta[MAXJ];
|
| 219 |
float genPhi[MAXJ];
|
| 220 |
float JECUnc[MAXJ];
|
| 221 |
float vtxMass[MAXJ];
|
| 222 |
float vtx3dL [MAXJ];
|
| 223 |
float vtx3deL[MAXJ];
|
| 224 |
} JetInfo;
|
| 225 |
|
| 226 |
int main(int argc, char* argv[])
|
| 227 |
{
|
| 228 |
gROOT->Reset();
|
| 229 |
|
| 230 |
TTree *_outTree;
|
| 231 |
METInfo MET;
|
| 232 |
MHTInfo MHT;
|
| 233 |
TopInfo top;
|
| 234 |
EventInfo EVENT;
|
| 235 |
// JetInfo jet1,jet2, addJet1, addJet2;
|
| 236 |
// lepton1,lepton2;
|
| 237 |
JetInfo hJets, aJets;
|
| 238 |
LeptonInfo vLeptons, aLeptons;
|
| 239 |
int naJets=0, nhJets=0;
|
| 240 |
TrackInfo H;
|
| 241 |
TrackInfo V;
|
| 242 |
int nvlep=0,nalep=0;
|
| 243 |
|
| 244 |
float jjdr,jjdPhi,jjdEta,HVdPhi,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMET, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight;
|
| 245 |
float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart;
|
| 246 |
|
| 247 |
int Vtype,numJets,numBJets,eventFlav;
|
| 248 |
// bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
|
| 249 |
bool triggerFlags[500],hbhe;
|
| 250 |
float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.;
|
| 251 |
// ----------------------------------------------------------------------
|
| 252 |
// First Part:
|
| 253 |
//
|
| 254 |
// * enable the AutoLibraryLoader
|
| 255 |
// * book the histograms of interest
|
| 256 |
// * open the input file
|
| 257 |
// ----------------------------------------------------------------------
|
| 258 |
|
| 259 |
// load framework libraries
|
| 260 |
gSystem->Load( "libFWCoreFWLite" );
|
| 261 |
gSystem->Load("libDataFormatsFWLite");
|
| 262 |
AutoLibraryLoader::enable();
|
| 263 |
|
| 264 |
// parse arguments
|
| 265 |
if ( argc < 2 ) {
|
| 266 |
return 0;
|
| 267 |
}
|
| 268 |
|
| 269 |
std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
|
| 270 |
std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
|
| 271 |
|
| 272 |
// get the python configuration
|
| 273 |
PythonProcessDesc builder(argv[1]);
|
| 274 |
const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
|
| 275 |
const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
|
| 276 |
const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
|
| 277 |
std::vector<edm::LuminosityBlockRange> jsonVector;
|
| 278 |
if ( in.exists("lumisToProcess") )
|
| 279 |
{
|
| 280 |
std::vector<edm::LuminosityBlockRange> const & lumisTemp =
|
| 281 |
in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
|
| 282 |
jsonVector.resize( lumisTemp.size() );
|
| 283 |
copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
|
| 284 |
}
|
| 285 |
|
| 286 |
// now get each parameter
|
| 287 |
int maxEvents_( in.getParameter<int>("maxEvents") );
|
| 288 |
unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
|
| 289 |
std::string outputFile_( out.getParameter<std::string>("fileName" ) );
|
| 290 |
std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
|
| 291 |
double btagThr = ana.getParameter<double>("bJetCountThreshold" );
|
| 292 |
bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
|
| 293 |
bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
|
| 294 |
bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
|
| 295 |
HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ );
|
| 296 |
HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW);
|
| 297 |
|
| 298 |
TriggerWeight triggerWeight(ana);
|
| 299 |
BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
|
| 300 |
BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
|
| 301 |
|
| 302 |
std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
|
| 303 |
// std::string inputFile( in.getParameter<std::string> ("fileName") );
|
| 304 |
|
| 305 |
|
| 306 |
std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
|
| 307 |
std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
|
| 308 |
|
| 309 |
bool isMC_( ana.getParameter<bool>("isMC") );
|
| 310 |
TriggerReader trigger(isMC_);
|
| 311 |
TriggerReader patFilters(false);
|
| 312 |
|
| 313 |
edm::LumiReWeighting lumiWeights;
|
| 314 |
if(isMC_)
|
| 315 |
{
|
| 316 |
lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
|
| 317 |
}
|
| 318 |
|
| 319 |
// TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
|
| 320 |
// TH1F * pu = new TH1F("pileup","",-0.5,24.5,25);
|
| 321 |
TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
|
| 322 |
TH1F * count = new TH1F("Count","Count", 1,0,2 );
|
| 323 |
TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
|
| 324 |
_outTree = new TTree("tree", "myTree");
|
| 325 |
|
| 326 |
_outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F");
|
| 327 |
_outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
|
| 328 |
_outTree->Branch("nhJets" , &nhJets , "nhJets/I");
|
| 329 |
_outTree->Branch("naJets" , &naJets , "naJets/I");
|
| 330 |
|
| 331 |
_outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
|
| 332 |
_outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
|
| 333 |
_outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
|
| 334 |
_outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
|
| 335 |
_outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
|
| 336 |
_outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
|
| 337 |
_outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
|
| 338 |
_outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
|
| 339 |
_outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
|
| 340 |
_outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
|
| 341 |
_outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
|
| 342 |
_outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
|
| 343 |
_outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
|
| 344 |
_outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
|
| 345 |
_outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
|
| 346 |
_outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
|
| 347 |
_outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
|
| 348 |
_outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[naJets]/F");
|
| 349 |
_outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[naJets]/F");
|
| 350 |
_outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[naJets]/F");
|
| 351 |
|
| 352 |
_outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
|
| 353 |
_outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
|
| 354 |
_outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
|
| 355 |
_outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
|
| 356 |
_outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
|
| 357 |
_outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
|
| 358 |
_outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
|
| 359 |
_outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
|
| 360 |
_outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
|
| 361 |
_outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
|
| 362 |
_outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
|
| 363 |
_outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
|
| 364 |
_outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
|
| 365 |
_outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
|
| 366 |
_outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
|
| 367 |
_outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
|
| 368 |
_outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
|
| 369 |
_outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
|
| 370 |
_outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
|
| 371 |
_outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
|
| 372 |
|
| 373 |
_outTree->Branch("jjdr" , &jjdr , "jjdr/F" );
|
| 374 |
_outTree->Branch("jjdPhi" , &jjdPhi , "jjdPhi/F" );
|
| 375 |
_outTree->Branch("jjdEta" , &jjdEta , "jjdEta/F" );
|
| 376 |
_outTree->Branch("numJets" , &numJets , "numJets/I" );
|
| 377 |
_outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
|
| 378 |
_outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
|
| 379 |
_outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
|
| 380 |
_outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
|
| 381 |
_outTree->Branch("genZpt" , &genZpt , "genZpt/F");
|
| 382 |
_outTree->Branch("genWpt" , &genWpt , "genWpt/F");
|
| 383 |
_outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
|
| 384 |
_outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
|
| 385 |
_outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
|
| 386 |
_outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
|
| 387 |
_outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
|
| 388 |
|
| 389 |
|
| 390 |
_outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
|
| 391 |
_outTree->Branch("PUweight", &PUweight , "PUweight/F");
|
| 392 |
_outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
|
| 393 |
|
| 394 |
|
| 395 |
|
| 396 |
|
| 397 |
_outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
|
| 398 |
_outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
|
| 399 |
_outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
|
| 400 |
_outTree->Branch("VMt" , &VMt , "VMt/F" );
|
| 401 |
|
| 402 |
_outTree->Branch("nvlep" , &nvlep , "nvlep/I");
|
| 403 |
_outTree->Branch("nalep" , &nalep , "nalep/I");
|
| 404 |
|
| 405 |
_outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
|
| 406 |
_outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
|
| 407 |
_outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
|
| 408 |
_outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
|
| 409 |
_outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
|
| 410 |
_outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
|
| 411 |
_outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
|
| 412 |
_outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
|
| 413 |
_outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
|
| 414 |
_outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
|
| 415 |
_outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
|
| 416 |
_outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
|
| 417 |
_outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
|
| 418 |
_outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
|
| 419 |
_outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
|
| 420 |
|
| 421 |
_outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
|
| 422 |
_outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
|
| 423 |
_outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
|
| 424 |
_outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
|
| 425 |
_outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
|
| 426 |
_outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
|
| 427 |
_outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
|
| 428 |
_outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
|
| 429 |
_outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
|
| 430 |
_outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
|
| 431 |
_outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
|
| 432 |
_outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
|
| 433 |
_outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
|
| 434 |
_outTree->Branch("aLepton_id80",aLeptons.id80 ,"id[nalep]/F");
|
| 435 |
_outTree->Branch("aLepton_id95",aLeptons.id95 ,"id[nalep]/F");
|
| 436 |
|
| 437 |
_outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
|
| 438 |
|
| 439 |
_outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
|
| 440 |
_outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
|
| 441 |
_outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
|
| 442 |
_outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
|
| 443 |
|
| 444 |
std::stringstream s;
|
| 445 |
s << "triggerFlags[" << triggers.size() << "]/b";
|
| 446 |
_outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
|
| 447 |
|
| 448 |
|
| 449 |
_outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
|
| 450 |
_outTree->Branch("hbhe" , &hbhe , "hbhe/b");
|
| 451 |
_outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
|
| 452 |
_outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
|
| 453 |
_outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
|
| 454 |
|
| 455 |
int ievt=0;
|
| 456 |
int totalcount=0;
|
| 457 |
|
| 458 |
// TFile* inFile = new TFile(inputFile.c_str(), "read");
|
| 459 |
for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
|
| 460 |
std::cout << iFile << std::endl;
|
| 461 |
TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
|
| 462 |
if(inFile==0) continue;
|
| 463 |
|
| 464 |
// loop the events
|
| 465 |
|
| 466 |
fwlite::Event ev(inFile);
|
| 467 |
for(ev.toBegin(); !ev.atEnd(); ++ev, ++ievt)
|
| 468 |
{
|
| 469 |
count->Fill(1.);
|
| 470 |
|
| 471 |
fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
|
| 472 |
vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 473 |
const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
|
| 474 |
PUweight=1.;
|
| 475 |
if(isMC_){
|
| 476 |
// PU weights
|
| 477 |
std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
|
| 478 |
PUweight = lumiWeights.weight( puit->second /3. );
|
| 479 |
}
|
| 480 |
countWithPU->Fill(PUweight);
|
| 481 |
|
| 482 |
//Write event info
|
| 483 |
EVENT.run = ev.id().run();
|
| 484 |
EVENT.lumi = ev.id().luminosityBlock();
|
| 485 |
EVENT.event = ev.id().event();
|
| 486 |
EVENT.json = jsonContainsEvent (jsonVector, ev);
|
| 487 |
|
| 488 |
|
| 489 |
const std::vector<VHbbCandidate> * candZ ;
|
| 490 |
const std::vector<VHbbCandidate> * candW ;
|
| 491 |
const VHbbEvent * iEvent =0;
|
| 492 |
if(fromCandidate)
|
| 493 |
{
|
| 494 |
fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
|
| 495 |
vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
|
| 496 |
candZ = vhbbCandHandleZ.product();
|
| 497 |
|
| 498 |
fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
|
| 499 |
vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
|
| 500 |
candW = vhbbCandHandle.product();
|
| 501 |
}
|
| 502 |
else
|
| 503 |
{
|
| 504 |
candZlocal->clear();
|
| 505 |
candWlocal->clear();
|
| 506 |
fwlite::Handle< VHbbEvent > vhbbHandle;
|
| 507 |
vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 508 |
iEvent = vhbbHandle.product();
|
| 509 |
algoZ->run(vhbbHandle.product(),*candZlocal);
|
| 510 |
algoW->run(vhbbHandle.product(),*candWlocal);
|
| 511 |
candZ= candZlocal;
|
| 512 |
candW= candWlocal;
|
| 513 |
|
| 514 |
/* for(size_t m=0;m<iEvent->muInfo.size();m++)
|
| 515 |
{
|
| 516 |
|
| 517 |
if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
|
| 518 |
fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
|
| 519 |
{
|
| 520 |
std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
|
| 521 |
}
|
| 522 |
}
|
| 523 |
*/
|
| 524 |
|
| 525 |
}
|
| 526 |
|
| 527 |
const std::vector<VHbbCandidate> * cand = candZ;
|
| 528 |
|
| 529 |
|
| 530 |
|
| 531 |
/* fwlite::Handle< VHbbEvent > vhbbHandle;
|
| 532 |
vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
|
| 533 |
const VHbbEvent iEvent = *vhbbHandle.product();
|
| 534 |
*/
|
| 535 |
|
| 536 |
// std::clog << "Filling tree "<< std::endl;
|
| 537 |
bool isW=false;
|
| 538 |
|
| 539 |
|
| 540 |
if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
|
| 541 |
if(cand->size() > 1 )
|
| 542 |
{
|
| 543 |
std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
|
| 544 |
}
|
| 545 |
if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
|
| 546 |
if(cand->size() == 0)
|
| 547 |
{
|
| 548 |
// std::cout << "W event loss due to tigther cuts" << std::endl;
|
| 549 |
continue;
|
| 550 |
}
|
| 551 |
const VHbbCandidate & vhCand = cand->at(0);
|
| 552 |
patFilters.setEvent(&ev,"VH");
|
| 553 |
hbhe = patFilters.accept("hbhe");
|
| 554 |
|
| 555 |
trigger.setEvent(&ev);
|
| 556 |
for(size_t j=0;j < triggers.size();j++)
|
| 557 |
triggerFlags[j]=trigger.accept(triggers[j]);
|
| 558 |
|
| 559 |
eventFlav=0;
|
| 560 |
if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
|
| 561 |
else if(aux.mcC.size() > 0) eventFlav=4;
|
| 562 |
|
| 563 |
|
| 564 |
H.mass = vhCand.H.p4.M();
|
| 565 |
H.pt = vhCand.H.p4.Pt();
|
| 566 |
H.eta = vhCand.H.p4.Eta();
|
| 567 |
H.phi = vhCand.H.p4.Phi();
|
| 568 |
V.mass = vhCand.V.p4.M();
|
| 569 |
V.pt = vhCand.V.p4.Pt();
|
| 570 |
V.eta = vhCand.V.p4.Eta();
|
| 571 |
V.phi = vhCand.V.p4.Phi();
|
| 572 |
nhJets=2;
|
| 573 |
hJets.set(vhCand.H.jets[0],0);
|
| 574 |
hJets.set(vhCand.H.jets[1],1);
|
| 575 |
aJets.reset();
|
| 576 |
naJets=vhCand.additionalJets.size();
|
| 577 |
numBJets=0;
|
| 578 |
if(vhCand.H.jets[0].csv> btagThr) numBJets++;
|
| 579 |
if(vhCand.H.jets[1].csv> btagThr) numBJets++;
|
| 580 |
for( int j=0; j < naJets && j < MAXJ; j++ )
|
| 581 |
{
|
| 582 |
aJets.set(vhCand.additionalJets[j],j);
|
| 583 |
if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
|
| 584 |
}
|
| 585 |
numJets = vhCand.additionalJets.size()+2;
|
| 586 |
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());
|
| 587 |
jjdPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
|
| 588 |
jjdEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
|
| 589 |
HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
|
| 590 |
HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
|
| 591 |
VMt = vhCand.Mt() ;
|
| 592 |
deltaPullAngle = vhCand.H.deltaTheta;
|
| 593 |
|
| 594 |
|
| 595 |
hJets.cosTheta[0]= vhCand.H.helicities[0];
|
| 596 |
hJets.cosTheta[1]= vhCand.H.helicities[1];
|
| 597 |
|
| 598 |
MET.et = vhCand.V.mets.at(0).p4.Pt();
|
| 599 |
MET.phi = vhCand.V.mets.at(0).p4.Phi();
|
| 600 |
MET.sumet = vhCand.V.mets.at(0).sumEt;
|
| 601 |
MET.sig = vhCand.V.mets.at(0).metSig;
|
| 602 |
|
| 603 |
|
| 604 |
if(!fromCandidate) {
|
| 605 |
MHT.mht = iEvent->mht.p4.Pt();
|
| 606 |
MHT.phi = iEvent->mht.p4.Phi();
|
| 607 |
MHT.ht = iEvent->mht.sumEt;
|
| 608 |
MHT.sig = iEvent->mht.metSig;
|
| 609 |
|
| 610 |
}
|
| 611 |
|
| 612 |
|
| 613 |
Vtype = vhCand.candidateType;
|
| 614 |
|
| 615 |
//Loop on jets
|
| 616 |
double maxBtag=-99999;
|
| 617 |
minDeltaPhijetMET = 999;
|
| 618 |
TLorentzVector bJet;
|
| 619 |
std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
|
| 620 |
std::vector<float> jet30eta;
|
| 621 |
std::vector<float> jet30pt;
|
| 622 |
if(fromCandidate)
|
| 623 |
{
|
| 624 |
//Loop on Higgs Jets
|
| 625 |
for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
|
| 626 |
if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
|
| 627 |
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()) < minDeltaPhijetMET)
|
| 628 |
{
|
| 629 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi());
|
| 630 |
jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
|
| 631 |
}
|
| 632 |
btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
|
| 633 |
}
|
| 634 |
//Loop on Additional Jets
|
| 635 |
for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
|
| 636 |
if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
|
| 637 |
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()) < minDeltaPhijetMET)
|
| 638 |
{
|
| 639 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi());
|
| 640 |
jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
|
| 641 |
}
|
| 642 |
if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
|
| 643 |
{
|
| 644 |
if(vhCand.additionalJets[j].p4.Pt() > 20)
|
| 645 |
btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
|
| 646 |
}
|
| 647 |
}
|
| 648 |
}
|
| 649 |
else
|
| 650 |
{
|
| 651 |
for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
|
| 652 |
if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
|
| 653 |
if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()) < minDeltaPhijetMET)
|
| 654 |
{
|
| 655 |
minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi());
|
| 656 |
jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
|
| 657 |
}
|
| 658 |
if(iEvent->simpleJets2[j].p4.Pt() > 30)
|
| 659 |
{
|
| 660 |
jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
|
| 661 |
jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
|
| 662 |
}
|
| 663 |
|
| 664 |
//For events made with highest CSV, all jets in the event should be taken into account for "tagging" SF (anti tagging is a mess)
|
| 665 |
// 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
|
| 666 |
if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
|
| 667 |
{
|
| 668 |
if(iEvent->simpleJets2[j].p4.Pt() > 20)
|
| 669 |
btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
|
| 670 |
}
|
| 671 |
}
|
| 672 |
|
| 673 |
//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
|
| 674 |
// by a criteria (pt of the dijet) that is not btag SF dependent
|
| 675 |
if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
|
| 676 |
for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
|
| 677 |
}
|
| 678 |
|
| 679 |
}
|
| 680 |
vLeptons.reset();
|
| 681 |
weightTrig = 1.; // better to default to 1
|
| 682 |
TLorentzVector leptonForTop;
|
| 683 |
size_t firstAddMu=0;
|
| 684 |
size_t firstAddEle=0;
|
| 685 |
if(Vtype == VHbbCandidate::Zmumu ){
|
| 686 |
vLeptons.set(vhCand.V.muons[0],0,13);
|
| 687 |
vLeptons.set(vhCand.V.muons[1],1,13);
|
| 688 |
float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
|
| 689 |
float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
|
| 690 |
float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
|
| 691 |
float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
|
| 692 |
weightTrig = cweightID * cweightTrig;
|
| 693 |
nvlep=2;
|
| 694 |
firstAddMu=2;
|
| 695 |
}
|
| 696 |
if( Vtype == VHbbCandidate::Zee ){
|
| 697 |
vLeptons.set(vhCand.V.electrons[0],0,11);
|
| 698 |
vLeptons.set(vhCand.V.electrons[1],1,11);
|
| 699 |
nvlep=2;
|
| 700 |
firstAddEle=2;
|
| 701 |
std::vector<float> pt,eta;
|
| 702 |
pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
|
| 703 |
pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
|
| 704 |
weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
|
| 705 |
triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
|
| 706 |
weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
|
| 707 |
weightTrig = weightEleTrigElePart * weightEleRecoAndId;
|
| 708 |
|
| 709 |
|
| 710 |
|
| 711 |
}
|
| 712 |
if(Vtype == VHbbCandidate::Wmun ){
|
| 713 |
leptonForTop=vhCand.V.muons[0].p4;
|
| 714 |
vLeptons.set(vhCand.V.muons[0],0,13);
|
| 715 |
float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
|
| 716 |
float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
|
| 717 |
float cweightTrig = weightTrig1;
|
| 718 |
weightTrig = cweightID * cweightTrig;
|
| 719 |
nvlep=1;
|
| 720 |
firstAddMu=1;
|
| 721 |
}
|
| 722 |
if( Vtype == VHbbCandidate::Wen ){
|
| 723 |
leptonForTop=vhCand.V.electrons[0].p4;
|
| 724 |
vLeptons.set(vhCand.V.electrons[0],0,11);
|
| 725 |
nvlep=1;
|
| 726 |
firstAddEle=1;
|
| 727 |
float weightMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
|
| 728 |
float weightV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
|
| 729 |
weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
|
| 730 |
weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30eta,jet30pt)*triggerWeight.scalePFMHTEle(MET.et);
|
| 731 |
weightEleTrigElePart= weightV4; //this is for debugging only, checking only the V4 part
|
| 732 |
|
| 733 |
weightMay*=weightEleRecoAndId;
|
| 734 |
weightV4*=weightEleRecoAndId;
|
| 735 |
weightV4*=weightEleTrigJetMETPart;
|
| 736 |
weightTrig = weightMay * 0.187 + weightV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
|
| 737 |
}
|
| 738 |
if( Vtype == VHbbCandidate::Znn ){
|
| 739 |
nvlep=0;
|
| 740 |
float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
|
| 741 |
weightTrig = weightTrig1;
|
| 742 |
weightTrigMET = weightTrig1;
|
| 743 |
}
|
| 744 |
|
| 745 |
aLeptons.reset();
|
| 746 |
nalep=0;
|
| 747 |
if(fromCandidate)
|
| 748 |
{
|
| 749 |
for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
|
| 750 |
for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
|
| 751 |
}
|
| 752 |
else
|
| 753 |
{
|
| 754 |
for(size_t j=0;j< iEvent->muInfo.size();j++)
|
| 755 |
{
|
| 756 |
if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
|
| 757 |
aLeptons.set(iEvent->muInfo[j],nalep++,13);
|
| 758 |
}
|
| 759 |
for(size_t j=0;j< iEvent->eleInfo.size();j++)
|
| 760 |
{
|
| 761 |
if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
|
| 762 |
aLeptons.set(iEvent->eleInfo[j],nalep++,11);
|
| 763 |
}
|
| 764 |
|
| 765 |
}
|
| 766 |
|
| 767 |
|
| 768 |
if(isMC_)
|
| 769 |
{
|
| 770 |
//std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
|
| 771 |
if ( btagJetInfos.size()< 10)
|
| 772 |
{
|
| 773 |
btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
|
| 774 |
btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
|
| 775 |
btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
|
| 776 |
}
|
| 777 |
else
|
| 778 |
{
|
| 779 |
std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
|
| 780 |
//TODO: revert to random throw for this cases
|
| 781 |
btag1T2CSF = 1.;
|
| 782 |
btag2TSF = 1.;
|
| 783 |
btag1TSF = 1.;
|
| 784 |
}
|
| 785 |
}
|
| 786 |
|
| 787 |
if(maxBtag > -99999)
|
| 788 |
{
|
| 789 |
TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
|
| 790 |
top.mass = topQuark.p4.M();
|
| 791 |
top.pt = topQuark.p4.Pt();
|
| 792 |
top.wMass = topQuark.p4W.M();
|
| 793 |
} else {
|
| 794 |
top.mass = -99;
|
| 795 |
top.pt = -99;
|
| 796 |
top.wMass = -99;
|
| 797 |
}
|
| 798 |
|
| 799 |
|
| 800 |
|
| 801 |
//FIXME: too much warnings... figure out why
|
| 802 |
// gendrcc=aux.genCCDeltaR();
|
| 803 |
// gendrbb=aux.genBBDeltaR();
|
| 804 |
genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
|
| 805 |
genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
|
| 806 |
|
| 807 |
|
| 808 |
/// Compute pull angle from AK7
|
| 809 |
if(!fromCandidate){
|
| 810 |
std::vector<VHbbEvent::SimpleJet> ak7 = iEvent->simpleJets3;
|
| 811 |
if(ak7.size() > 1){
|
| 812 |
CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
|
| 813 |
CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
|
| 814 |
std::vector<VHbbEvent::SimpleJet> ak7_matched;
|
| 815 |
|
| 816 |
std::sort( ak7.begin(),ak7.end(),deltaRComparatorJ1 );
|
| 817 |
ak7_matched.push_back(ak7[0]);
|
| 818 |
|
| 819 |
if(ak7[1].p4.DeltaR(vhCand.H.jets[0].p4) > 0.5)
|
| 820 |
ak7.erase(ak7.begin());
|
| 821 |
|
| 822 |
std::sort( ak7.begin(),ak7.end(),deltaRComparatorJ2 );
|
| 823 |
if( abs(ak7[0].p4.Pt() - ak7_matched[1].p4.Pt()) > 0.1 )
|
| 824 |
ak7_matched.push_back(ak7[0]);
|
| 825 |
else
|
| 826 |
ak7_matched.push_back(ak7[1]);
|
| 827 |
|
| 828 |
CompareJetPt ptComparator;
|
| 829 |
std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
|
| 830 |
if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
|
| 831 |
and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5 ){
|
| 832 |
deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
|
| 833 |
}
|
| 834 |
}
|
| 835 |
}
|
| 836 |
|
| 837 |
|
| 838 |
|
| 839 |
|
| 840 |
|
| 841 |
_outTree->Fill();
|
| 842 |
|
| 843 |
}// closed event loop
|
| 844 |
|
| 845 |
std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
|
| 846 |
inFile->Close();
|
| 847 |
// close input file
|
| 848 |
} // loop on files
|
| 849 |
|
| 850 |
|
| 851 |
std::cout << "Events: " << ievt <<std::endl;
|
| 852 |
std::cout << "TotalCount: " << totalcount <<std::endl;
|
| 853 |
|
| 854 |
|
| 855 |
|
| 856 |
_outFile->cd();
|
| 857 |
|
| 858 |
_outTree->Write();
|
| 859 |
_outFile->Write();
|
| 860 |
_outFile->Close();
|
| 861 |
return 0;
|
| 862 |
}
|
| 863 |
|
| 864 |
|
