| 1 |
/* A macro for making histograms for an LJMET analysis
|
| 2 |
*/
|
| 3 |
|
| 4 |
#include "TFile.h"
|
| 5 |
#include "TTree.h"
|
| 6 |
#include "TBranch.h"
|
| 7 |
#include "TH1.h"
|
| 8 |
#include "TH2.h"
|
| 9 |
#include "TCanvas.h"
|
| 10 |
#include "TLegend.h"
|
| 11 |
#include "TSystem.h"
|
| 12 |
#include "TLorentzVector.h"
|
| 13 |
#include "TVector.h"
|
| 14 |
#include "TMath.h"
|
| 15 |
|
| 16 |
#include "PhysicsTools/SelectorUtils/interface/strbitset.h"
|
| 17 |
#include "DataFormats/PatCandidates/interface/Jet.h"
|
| 18 |
#include "DataFormats/PatCandidates/interface/MET.h"
|
| 19 |
#include "DataFormats/PatCandidates/interface/Photon.h"
|
| 20 |
#include "DataFormats/PatCandidates/interface/Muon.h"
|
| 21 |
#include "DataFormats/PatCandidates/interface/Electron.h"
|
| 22 |
#include "DataFormats/PatCandidates/interface/Tau.h"
|
| 23 |
#include "PhysicsTools/SelectorUtils/interface/Selector.h"
|
| 24 |
#include "LJMet/Com/interface/MuonSelectionFunctor.h"
|
| 25 |
#include "LJMet/Com/interface/JetIDSelectionFunctor.h"
|
| 26 |
#include "LJMet/Com/interface/MetSelectionFunctor.h"
|
| 27 |
#include "FWCore/FWLite/interface/AutoLibraryLoader.h"
|
| 28 |
#include "PhysicsTools/FWLite/interface/TFileService.h"
|
| 29 |
#include "DataFormats/FWLite/interface/ChainEvent.h"
|
| 30 |
#include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
|
| 31 |
#include "FWCore/ParameterSet/interface/ProcessDesc.h"
|
| 32 |
#include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
|
| 33 |
#include "DataFormats/VertexReco/interface/Vertex.h"
|
| 34 |
#include "DataFormats/TrackReco/interface/Track.h"
|
| 35 |
#include "DataFormats/Candidate/interface/LeafCandidate.h"
|
| 36 |
#include "DataFormats/TrackReco/interface/TrackBase.h"
|
| 37 |
|
| 38 |
|
| 39 |
|
| 40 |
#include <iostream>
|
| 41 |
#include <cmath> //necessary for absolute function fabs()
|
| 42 |
|
| 43 |
using namespace std;
|
| 44 |
|
| 45 |
///-------------------------
|
| 46 |
/// DRIVER FUNCTION
|
| 47 |
///-------------------------
|
| 48 |
|
| 49 |
// -*- C++ -*-
|
| 50 |
|
| 51 |
// CMS includes
|
| 52 |
#include "FWCore/Utilities/interface/InputTag.h"
|
| 53 |
#include "DataFormats/Common/interface/Handle.h"
|
| 54 |
#include "DataFormats/Common/interface/Ptr.h"
|
| 55 |
#include "DataFormats/PatCandidates/interface/Jet.h"
|
| 56 |
#include "PhysicsTools/SelectorUtils/interface/EventSelector.h"
|
| 57 |
|
| 58 |
|
| 59 |
// Root includes
|
| 60 |
#include "TROOT.h"
|
| 61 |
|
| 62 |
using namespace std;
|
| 63 |
|
| 64 |
|
| 65 |
class LJMetStudiesSelector : public EventSelector {
|
| 66 |
public:
|
| 67 |
LJMetStudiesSelector( edm::ParameterSet const & muonSelParams,
|
| 68 |
edm::ParameterSet const & jetSelParams,
|
| 69 |
edm::ParameterSet const & metSelParams,
|
| 70 |
edm::ParameterSet const & params ) :
|
| 71 |
muonSel_ (new MuonSelectionFunctor(muonSelParams)),
|
| 72 |
jetSel_ (new JetIDSelectionFunctor(jetSelParams)),
|
| 73 |
metSel_ (new MetSelectionFunctor(metSelParams)),
|
| 74 |
muonSrc_ (params.getParameter<edm::InputTag>("muonSrc")),
|
| 75 |
jetSrc_ (params.getParameter<edm::InputTag>("jetSrc")),
|
| 76 |
metSrc_ (params.getParameter<edm::InputTag>("metSrc"))
|
| 77 |
{
|
| 78 |
|
| 79 |
push_back("Jet cuts");
|
| 80 |
push_back("Min jet multiplicity", 2);
|
| 81 |
|
| 82 |
set("Jet cuts", true);
|
| 83 |
set("Min jet multiplicity", 2);
|
| 84 |
|
| 85 |
push_back("Muon Cuts");
|
| 86 |
push_back("Min muon multiplicity", 1);
|
| 87 |
|
| 88 |
set("Muon Cuts", true);
|
| 89 |
set("Min muon multiplicity", 1);
|
| 90 |
|
| 91 |
|
| 92 |
}
|
| 93 |
|
| 94 |
virtual ~LJMetStudiesSelector() {}
|
| 95 |
|
| 96 |
virtual bool operator()( edm::EventBase const & event, std::strbitset & ret){
|
| 97 |
|
| 98 |
//cout << "LJMetStudiesSelector()" << endl;
|
| 99 |
|
| 100 |
std::strbitset retMuon = muonSel_->getBitTemplate();
|
| 101 |
std::strbitset retJet = jetSel_->getBitTemplate();
|
| 102 |
std::strbitset retMet = metSel_->getBitTemplate();
|
| 103 |
|
| 104 |
while(1){
|
| 105 |
|
| 106 |
if ( considerCut("Muon Cuts") ) {
|
| 107 |
passCut(ret, "Muon Cuts");
|
| 108 |
event.getByLabel( muonSrc_, h_muons_ );
|
| 109 |
}
|
| 110 |
else break;
|
| 111 |
|
| 112 |
int nMuons = 0;
|
| 113 |
for ( vector<pat::Muon>::const_iterator mu = h_muons_->begin();
|
| 114 |
mu != h_muons_->end(); mu++){
|
| 115 |
|
| 116 |
retMuon.set(false);
|
| 117 |
bool pass = (*muonSel_)( *mu, retMuon );
|
| 118 |
|
| 119 |
if ( pass ){
|
| 120 |
++nMuons;
|
| 121 |
|
| 122 |
// save the first and second good muon
|
| 123 |
if (nMuons == 1) muon0_ = edm::Ptr<pat::Muon>( h_muons_, 0);
|
| 124 |
if (nMuons == 2) muon1_ = edm::Ptr<pat::Muon>( h_muons_, 1);
|
| 125 |
if (nMuons == 3) muon2_ = edm::Ptr<pat::Muon>( h_muons_, 2);
|
| 126 |
if (nMuons == 4) muon3_ = edm::Ptr<pat::Muon>( h_muons_, 3);
|
| 127 |
if (nMuons == 5) muon4_ = edm::Ptr<pat::Muon>( h_muons_, 4);
|
| 128 |
if (nMuons == 6) muon5_ = edm::Ptr<pat::Muon>( h_muons_, 5);
|
| 129 |
|
| 130 |
|
| 131 |
}
|
| 132 |
}
|
| 133 |
|
| 134 |
if( nMuons >= cut("Min muon multiplicity", int()) ) passCut(ret, "Min muon multiplicity");
|
| 135 |
else break;
|
| 136 |
|
| 137 |
//
|
| 138 |
//_____ Jet cuts __________________________________
|
| 139 |
//
|
| 140 |
|
| 141 |
|
| 142 |
if ( considerCut("Jet cuts") ) {
|
| 143 |
passCut(ret, "Jet cuts");
|
| 144 |
event.getByLabel( metSrc_, h_met_ );
|
| 145 |
event.getByLabel( jetSrc_, h_jets_ );
|
| 146 |
event.getByLabel( muonSrc_, h_muons_ );
|
| 147 |
}
|
| 148 |
else break;
|
| 149 |
|
| 150 |
// loop over jets
|
| 151 |
int nJets = 0;
|
| 152 |
for (vector<pat::Jet>::const_iterator jet = h_jets_->begin();
|
| 153 |
jet != h_jets_->end(); jet++){
|
| 154 |
|
| 155 |
retJet.set(false);
|
| 156 |
bool pass = (*jetSel_)( *jet, retJet );
|
| 157 |
|
| 158 |
if ( pass ){
|
| 159 |
++nJets;
|
| 160 |
|
| 161 |
// save the first two good jets
|
| 162 |
if (nJets == 1) jet0_ = edm::Ptr<pat::Jet>( h_jets_, 0);
|
| 163 |
if (nJets == 2) jet1_ = edm::Ptr<pat::Jet>( h_jets_, 1);
|
| 164 |
if (nJets == 3) jet2_ = edm::Ptr<pat::Jet>( h_jets_, 2);
|
| 165 |
if (nJets == 4) jet3_ = edm::Ptr<pat::Jet>( h_jets_, 3);
|
| 166 |
if (nJets == 5) jet4_ = edm::Ptr<pat::Jet>( h_jets_, 4);
|
| 167 |
if (nJets == 6) jet5_ = edm::Ptr<pat::Jet>( h_jets_, 5);
|
| 168 |
if (nJets == 7) jet6_ = edm::Ptr<pat::Jet>( h_jets_, 6);
|
| 169 |
|
| 170 |
}
|
| 171 |
}
|
| 172 |
if ( nJets >= cut("Min jet multiplicity",int()) ) passCut(ret, "Min jet multiplicity");
|
| 173 |
else break;
|
| 174 |
|
| 175 |
break;
|
| 176 |
|
| 177 |
} //end of while loop
|
| 178 |
|
| 179 |
|
| 180 |
|
| 181 |
return (bool)ret;
|
| 182 |
|
| 183 |
setIgnored(ret);
|
| 184 |
|
| 185 |
return false;
|
| 186 |
}// end of method
|
| 187 |
|
| 188 |
|
| 189 |
|
| 190 |
|
| 191 |
|
| 192 |
boost::shared_ptr<MetSelectionFunctor> const & metSel() const { return metSel_;}
|
| 193 |
// FIXME: implement proper selection functors
|
| 194 |
boost::shared_ptr<JetIDSelectionFunctor> const & jetSel() const { return jetSel_;}
|
| 195 |
boost::shared_ptr<MuonSelectionFunctor> const & muonSel() const { return muonSel_;}
|
| 196 |
|
| 197 |
|
| 198 |
pat::MET const & met() const { return *met_; }
|
| 199 |
|
| 200 |
vector<pat::Jet> const & allJets () const { return *h_jets_; }
|
| 201 |
pat::Jet const & jet0() const { return *jet0_; }
|
| 202 |
pat::Jet const & jet1() const { return *jet1_; }
|
| 203 |
pat::Jet const & jet2() const { return *jet2_; }
|
| 204 |
pat::Jet const & jet3() const { return *jet3_; }
|
| 205 |
pat::Jet const & jet4() const { return *jet4_; }
|
| 206 |
pat::Jet const & jet5() const { return *jet5_; }
|
| 207 |
pat::Jet const & jet6() const { return *jet6_; }
|
| 208 |
|
| 209 |
|
| 210 |
vector<pat::Muon> const & allMuons () const { return *h_muons_; }
|
| 211 |
pat::Muon const & muon0() const { return *muon0_; }
|
| 212 |
pat::Muon const & muon1() const { return *muon1_; }
|
| 213 |
pat::Muon const & muon2() const { return *muon2_; }
|
| 214 |
pat::Muon const & muon3() const { return *muon3_; }
|
| 215 |
pat::Muon const & muon4() const { return *muon4_; }
|
| 216 |
pat::Muon const & muon5() const { return *muon5_; }
|
| 217 |
|
| 218 |
|
| 219 |
|
| 220 |
protected:
|
| 221 |
boost::shared_ptr<MuonSelectionFunctor> muonSel_;
|
| 222 |
boost::shared_ptr<JetIDSelectionFunctor> jetSel_;
|
| 223 |
boost::shared_ptr<MetSelectionFunctor> metSel_;
|
| 224 |
|
| 225 |
edm::InputTag muonSrc_;
|
| 226 |
edm::InputTag jetSrc_;
|
| 227 |
edm::InputTag metSrc_;
|
| 228 |
|
| 229 |
edm::Handle<vector<pat::MET> > h_met_;
|
| 230 |
edm::Handle<vector<pat::Jet> > h_jets_;
|
| 231 |
edm::Handle<vector<pat::Muon> > h_muons_;
|
| 232 |
|
| 233 |
edm::Ptr<pat::MET> met_;
|
| 234 |
edm::Ptr<pat::Jet> jet0_;
|
| 235 |
edm::Ptr<pat::Jet> jet1_;
|
| 236 |
edm::Ptr<pat::Jet> jet2_;
|
| 237 |
edm::Ptr<pat::Jet> jet3_;
|
| 238 |
edm::Ptr<pat::Jet> jet4_;
|
| 239 |
edm::Ptr<pat::Jet> jet5_;
|
| 240 |
edm::Ptr<pat::Jet> jet6_;
|
| 241 |
edm::Ptr<pat::Muon> muon0_;
|
| 242 |
edm::Ptr<pat::Muon> muon1_;
|
| 243 |
edm::Ptr<pat::Muon> muon2_;
|
| 244 |
edm::Ptr<pat::Muon> muon3_;
|
| 245 |
edm::Ptr<pat::Muon> muon4_;
|
| 246 |
edm::Ptr<pat::Muon> muon5_;
|
| 247 |
};
|
| 248 |
|
| 249 |
///////////////////////////
|
| 250 |
// ///////////////////// //
|
| 251 |
// // Main Subroutine // //
|
| 252 |
// ///////////////////// //
|
| 253 |
///////////////////////////
|
| 254 |
|
| 255 |
int main (int argc, char* argv[])
|
| 256 |
{
|
| 257 |
|
| 258 |
|
| 259 |
if ( argc < 2 ) {
|
| 260 |
std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
|
| 261 |
return 0;
|
| 262 |
}
|
| 263 |
|
| 264 |
cout << "Hello from " << argv[0] << "!" << endl;
|
| 265 |
|
| 266 |
|
| 267 |
// load framework libraries
|
| 268 |
gSystem->Load( "libFWCoreFWLite" );
|
| 269 |
AutoLibraryLoader::enable();
|
| 270 |
|
| 271 |
|
| 272 |
cout << "Getting parameters" << endl;
|
| 273 |
// Get the python configuration
|
| 274 |
PythonProcessDesc builder(argv[1]);
|
| 275 |
boost::shared_ptr<edm::ProcessDesc> b = builder.processDesc();
|
| 276 |
boost::shared_ptr<edm::ParameterSet> parameters = b->getProcessPSet();
|
| 277 |
parameters->registerIt();
|
| 278 |
|
| 279 |
edm::ParameterSet const& ljmetStudiesParams = parameters->getParameter<edm::ParameterSet>("ljmetStudies");
|
| 280 |
edm::ParameterSet const& muonSelParameters = parameters->getParameter<edm::ParameterSet>("muonSelector");
|
| 281 |
edm::ParameterSet const& jetSelParameters = parameters->getParameter<edm::ParameterSet>("jetIDSelector");
|
| 282 |
edm::ParameterSet const& metSelParameters = parameters->getParameter<edm::ParameterSet>("metSelector");
|
| 283 |
edm::ParameterSet const& plotParameters = parameters->getParameter<edm::ParameterSet>("plotParameters");
|
| 284 |
edm::ParameterSet const& inputs = parameters->getParameter<edm::ParameterSet>("inputs");
|
| 285 |
edm::ParameterSet const& outputs = parameters->getParameter<edm::ParameterSet>("outputs");
|
| 286 |
|
| 287 |
cout << "setting up TFileService" << endl;
|
| 288 |
// book a set of histograms
|
| 289 |
fwlite::TFileService fs = fwlite::TFileService( outputs.getParameter<std::string>("outputName") );
|
| 290 |
TFileDirectory theDir = fs.mkdir( "histos" );
|
| 291 |
|
| 292 |
//TH1F * dummyhist = fs.make<TH1F>( "dummyHist", "Dummy histogram", 10, 0, 10 );
|
| 293 |
//dummyhist->Fill(1);
|
| 294 |
|
| 295 |
cout << "Setting up chain event" << endl;
|
| 296 |
// This object 'event' is used both to get all information from the
|
| 297 |
// event as well as to store histograms, etc.
|
| 298 |
fwlite::ChainEvent ev ( inputs.getParameter<std::vector<std::string> > ("fileNames") );
|
| 299 |
|
| 300 |
|
| 301 |
cout << "Booking histograms" << endl;
|
| 302 |
// Book histograms
|
| 303 |
|
| 304 |
std::map<std::string, TH1*> hists;
|
| 305 |
std::map<std::string, TTree*> tree;
|
| 306 |
|
| 307 |
|
| 308 |
int nMuons;
|
| 309 |
double px, py, pz, et_mu0, phi_mu0, pt0, phi0, pt_jet0, pt_jet1, pt_jet2, pt_jet3, pt_jet4 ;
|
| 310 |
double px_jet0, px_jet1, px_jet2, px_jet3, px_jet4 ;
|
| 311 |
double py_jet0, py_jet1, py_jet2, py_jet3, py_jet4 ;
|
| 312 |
double pz_jet0, pz_jet1, pz_jet2, pz_jet3, pz_jet4 ;
|
| 313 |
double eta_jet0, eta_jet1, eta_jet2, eta_jet3, eta_jet4 ;
|
| 314 |
double phi_jet0, phi_jet1, phi_jet2, phi_jet3, phi_jet4 ;
|
| 315 |
double pt_mu0, pt_mu1, pt_mu2, pt_mu3 ;
|
| 316 |
float bDiscriminator_tche_jet0;
|
| 317 |
float bDiscriminator_tchp_jet0;
|
| 318 |
float bDiscriminator_ssv_jet0;
|
| 319 |
float bDiscriminator_ssv_hp_jet0;
|
| 320 |
float bDiscriminator_tche_jet1;
|
| 321 |
float bDiscriminator_tchp_jet1;
|
| 322 |
float bDiscriminator_ssv_jet1;
|
| 323 |
float bDiscriminator_ssv_hp_jet1;
|
| 324 |
float bDiscriminator_tche_jet2;
|
| 325 |
float bDiscriminator_tchp_jet2;
|
| 326 |
float bDiscriminator_ssv_jet2;
|
| 327 |
float bDiscriminator_ssv_hp_jet2;
|
| 328 |
float bDiscriminator_tche_jet3;
|
| 329 |
float bDiscriminator_tchp_jet3;
|
| 330 |
float bDiscriminator_ssv_jet3;
|
| 331 |
float bDiscriminator_ssv_hp_jet3;
|
| 332 |
float bDiscriminator_tche_jet4;
|
| 333 |
float bDiscriminator_tchp_jet4;
|
| 334 |
float bDiscriminator_ssv_jet4;
|
| 335 |
float bDiscriminator_ssv_hp_jet4;
|
| 336 |
px = 0;
|
| 337 |
py = 0;
|
| 338 |
pz = 0;
|
| 339 |
|
| 340 |
bDiscriminator_tche_jet0 = 0;
|
| 341 |
bDiscriminator_tchp_jet0 = 0;
|
| 342 |
bDiscriminator_ssv_jet0 = 0;
|
| 343 |
bDiscriminator_ssv_hp_jet0 = 0;
|
| 344 |
bDiscriminator_tche_jet1 = 0;
|
| 345 |
bDiscriminator_tchp_jet1 = 0;
|
| 346 |
bDiscriminator_ssv_jet1 = 0;
|
| 347 |
bDiscriminator_ssv_hp_jet1 = 0;
|
| 348 |
bDiscriminator_tche_jet2 = 0;
|
| 349 |
bDiscriminator_tchp_jet2 = 0;
|
| 350 |
bDiscriminator_ssv_jet2 = 0;
|
| 351 |
bDiscriminator_ssv_hp_jet2 = 0;
|
| 352 |
bDiscriminator_tche_jet3 = 0;
|
| 353 |
bDiscriminator_tchp_jet3 = 0;
|
| 354 |
bDiscriminator_ssv_jet3 = 0;
|
| 355 |
bDiscriminator_ssv_hp_jet3 = 0;
|
| 356 |
bDiscriminator_tche_jet4 = 0;
|
| 357 |
bDiscriminator_tchp_jet4 = 0;
|
| 358 |
bDiscriminator_ssv_jet4 = 0;
|
| 359 |
bDiscriminator_ssv_hp_jet4 = 0;
|
| 360 |
|
| 361 |
tree["Kinematic_Variables" ] = fs.make<TTree>("Kinematic_Variables", "Kinematic_Variables", 64000000);
|
| 362 |
TBranch *b_nMuons = tree["Kinematic_Variables"]->Branch("nMuons", &nMuons, "nMuons/I");
|
| 363 |
TBranch *b_phi_Muon0 = tree["Kinematic_Variables"]->Branch("phi_mu0", &phi_mu0, "phi_mu0/D");
|
| 364 |
TBranch *b_pt_Muon0 = tree["Kinematic_Variables"]->Branch("pt_mu0", &pt_mu0, "pt_mu0/D");
|
| 365 |
TBranch *b_pt_Muon1 = tree["Kinematic_Variables"]->Branch("pt_mu1", &pt_mu1, "pt_mu1/D");
|
| 366 |
TBranch *b_pt_Muon2 = tree["Kinematic_Variables"]->Branch("pt_mu2", &pt_mu2, "pt_mu2/D");
|
| 367 |
TBranch *b_pt_Muon3 = tree["Kinematic_Variables"]->Branch("pt_mu3", &pt_mu3, "pt_mu3/D");
|
| 368 |
TBranch *b_pt_jet0 = tree["Kinematic_Variables"]->Branch("pt_jet0", &pt_jet0, "pt_jet0/D");
|
| 369 |
TBranch *b_pt_jet1 = tree["Kinematic_Variables"]->Branch("pt_jet1", &pt_jet1, "pt_jet1/D");
|
| 370 |
TBranch *b_pt_jet2 = tree["Kinematic_Variables"]->Branch("pt_jet2", &pt_jet2, "pt_jet2/D");
|
| 371 |
TBranch *b_pt_jet3 = tree["Kinematic_Variables"]->Branch("pt_jet3", &pt_jet3, "pt_jet3/D");
|
| 372 |
TBranch *b_pt_jet4 = tree["Kinematic_Variables"]->Branch("pt_jet4", &pt_jet4, "pt_jet4/D");
|
| 373 |
TBranch *b_px_jet0 = tree["Kinematic_Variables"]->Branch("px_jet0", &px_jet0, "px_jet0/D");
|
| 374 |
TBranch *b_px_jet1 = tree["Kinematic_Variables"]->Branch("px_jet1", &px_jet1, "px_jet1/D");
|
| 375 |
TBranch *b_px_jet2 = tree["Kinematic_Variables"]->Branch("px_jet2", &px_jet2, "px_jet2/D");
|
| 376 |
TBranch *b_px_jet3 = tree["Kinematic_Variables"]->Branch("px_jet3", &px_jet3, "px_jet3/D");
|
| 377 |
TBranch *b_px_jet4 = tree["Kinematic_Variables"]->Branch("px_jet4", &px_jet4, "px_jet4/D");
|
| 378 |
TBranch *b_py_jet0 = tree["Kinematic_Variables"]->Branch("py_jet0", &py_jet0, "py_jet0/D");
|
| 379 |
TBranch *b_py_jet1 = tree["Kinematic_Variables"]->Branch("py_jet1", &py_jet1, "py_jet1/D");
|
| 380 |
TBranch *b_py_jet2 = tree["Kinematic_Variables"]->Branch("py_jet2", &py_jet2, "py_jet2/D");
|
| 381 |
TBranch *b_py_jet3 = tree["Kinematic_Variables"]->Branch("py_jet3", &py_jet3, "py_jet3/D");
|
| 382 |
TBranch *b_py_jet4 = tree["Kinematic_Variables"]->Branch("py_jet4", &py_jet4, "py_jet4/D");
|
| 383 |
TBranch *b_pz_jet0 = tree["Kinematic_Variables"]->Branch("pz_jet0", &py_jet0, "pz_jet0/D");
|
| 384 |
TBranch *b_pz_jet1 = tree["Kinematic_Variables"]->Branch("pz_jet1", &py_jet1, "pz_jet1/D");
|
| 385 |
TBranch *b_pz_jet2 = tree["Kinematic_Variables"]->Branch("pz_jet2", &py_jet2, "pz_jet2/D");
|
| 386 |
TBranch *b_pz_jet3 = tree["Kinematic_Variables"]->Branch("pz_jet3", &py_jet3, "pz_jet3/D");
|
| 387 |
TBranch *b_pz_jet4 = tree["Kinematic_Variables"]->Branch("pz_jet4", &py_jet4, "pz_jet4/D");
|
| 388 |
TBranch *b_bDiscriminator_ssv_jet0 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_jet0", &bDiscriminator_ssv_jet0,
|
| 389 |
"bDiscriminator_ssv_jet0/F");
|
| 390 |
TBranch *b_bDiscriminator_ssv_hp_jet0 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_hp_jet0", &bDiscriminator_ssv_hp_jet0,
|
| 391 |
"bDiscriminator_ssv_hp_jet0/F");
|
| 392 |
TBranch *b_bDiscriminator_tche_jet0 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tche_jet0", &bDiscriminator_tche_jet0,"bDiscriminator_tche_jet0/F");
|
| 393 |
TBranch *b_bDiscriminator_tchp_jet0 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tchp_jet0", &bDiscriminator_tchp_jet0,"bDiscriminator_tchp_jet0/F");
|
| 394 |
TBranch *b_bDiscriminator_ssv_jet1 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_jet1", &bDiscriminator_ssv_jet1,
|
| 395 |
"bDiscriminator_ssv_jet1/F");
|
| 396 |
TBranch *b_bDiscriminator_ssv_hp_jet1 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_hp_jet1", &bDiscriminator_ssv_hp_jet1,
|
| 397 |
"bDiscriminator_ssv_hp_jet1/F");
|
| 398 |
TBranch *b_bDiscriminator_tche_jet1 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tche_jet1", &bDiscriminator_tche_jet1,"bDiscriminator_tche_jet1/F");
|
| 399 |
TBranch *b_bDiscriminator_tchp_jet1 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tchp_jet1", &bDiscriminator_tchp_jet1,"bDiscriminator_tche_jet1/F");
|
| 400 |
TBranch *b_bDiscriminator_ssv_jet2 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_jet2", &bDiscriminator_ssv_jet2,
|
| 401 |
"bDiscriminator_ssv_jet2/F");
|
| 402 |
TBranch *b_bDiscriminator_ssv_hp_jet2 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_hp_jet2", &bDiscriminator_ssv_hp_jet2,
|
| 403 |
"bDiscriminator_ssv_hp_jet2/F");
|
| 404 |
TBranch *b_bDiscriminator_tche_jet2 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tche_jet2", &bDiscriminator_tche_jet2,"bDiscriminator_tche_jet2/F");
|
| 405 |
TBranch *b_bDiscriminator_tchp_jet2 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tchp_jet2", &bDiscriminator_tchp_jet2,"bDiscriminator_tche_jet2/F");
|
| 406 |
TBranch *b_bDiscriminator_ssv_jet3 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_jet3", &bDiscriminator_ssv_jet3,
|
| 407 |
"bDiscriminator_ssv_jet3/F");
|
| 408 |
TBranch *b_bDiscriminator_ssv_hp_jet3 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_hp_jet3", &bDiscriminator_ssv_hp_jet3,
|
| 409 |
"bDiscriminator_ssv_hp_jet3/F");
|
| 410 |
TBranch *b_bDiscriminator_tche_jet3 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tche_jet3", &bDiscriminator_tche_jet3,"bDiscriminator_tche_jet3/F");
|
| 411 |
TBranch *b_bDiscriminator_tchp_jet3 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tchp_jet3", &bDiscriminator_tchp_jet3,"bDiscriminator_tche_jet3/F");
|
| 412 |
TBranch *b_bDiscriminator_ssv_jet4 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_jet4", &bDiscriminator_ssv_jet4,
|
| 413 |
"bDiscriminator_ssv_jet4/F");
|
| 414 |
TBranch *b_bDiscriminator_ssv_hp_jet4 = tree["Kinematic_Variables"]->Branch("bDiscriminator_ssv_hp_jet4", &bDiscriminator_ssv_hp_jet4,
|
| 415 |
"bDiscriminator_ssv_hp_jet4/F");
|
| 416 |
TBranch *b_bDiscriminator_tche_jet4 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tche_jet4", &bDiscriminator_tche_jet4,"bDiscriminator_tche_jet4/F");
|
| 417 |
TBranch *b_bDiscriminator_tchp_jet4 = tree["Kinematic_Variables"]->Branch("bDiscriminator_tchp_jet4", &bDiscriminator_tchp_jet4,"bDiscriminator_tche_jet4/F");
|
| 418 |
|
| 419 |
|
| 420 |
hists["hist_jet0Pt" ] = theDir.make<TH1D>( "hist_jet0Pt" , "Jet0 p_{T}", 400, 0, 400 ) ;
|
| 421 |
hists["hist_jet0Phi" ] = theDir.make<TH1D>( "hist_jet0Phi" , "MET phi", 400, -3, 3 ) ;
|
| 422 |
hists["hist_muon0Pt" ] = theDir.make<TH1D>( "hist_muon0Pt" , "Muon p_{T}", 400, 0, 400);
|
| 423 |
hists["hist_muon0Phi" ] = theDir.make<TH1D>( "hist_muon0Phi" , "Muon phi", 400, -3, 3);
|
| 424 |
hists["hist_jet1Pt" ] = theDir.make<TH1D>( "hist_jet1Pt" , "Jet1 p_{T}", 400, 0, 400);
|
| 425 |
hists["hist_jet0px" ] = theDir.make<TH1D>( "hist_jet0px" , "Jet0 px", 400, 0, 400);
|
| 426 |
hists["hist_jet0py" ] = theDir.make<TH1D>( "hist_jet0py" , "Jet0 py", 400, 0, 400);
|
| 427 |
hists["hist_jet1px" ] = theDir.make<TH1D>( "hist_jet1px" , "Jet1 px", 400, 0, 400);
|
| 428 |
hists["hist_jet1py" ] = theDir.make<TH1D>( "hist_jet1py" , "Jet1 py", 400, 0, 400);
|
| 429 |
hists["hist_cos_angle_jets01" ] = theDir.make<TH1D>( "hist_cos_angle_jets01" , "Cosine of the angle between jet 0 and 1", 40, -2, 2);
|
| 430 |
hists["hist_nvertices" ] = theDir.make<TH1D>( "hist_nvertices" , "Number of secondary vertices", 10, -0.5, 9.5);
|
| 431 |
hists["hist_b_discriminator_tche" ] = theDir.make<TH1D>( "hist_b_discriminator_tche" , "B Discriminator for Track Counting High Efficiency", 100, -50, 50);
|
| 432 |
hists["hist_b_discriminator_ssv" ] = theDir.make<TH1D>("hist_b_discriminator_ssv" , "B Discriminator for Simple Secondary Vertex Algorithm", 100, 0.0, 100);
|
| 433 |
hists["hist_b_discriminator_tchp" ] = theDir.make<TH1D>("hist_b_discriminator_tchp" , "B Discriminator for Track Counting High Purity", 100, -50, 50) ;
|
| 434 |
hists["hist_b_discriminator_ssv_hp" ] = theDir.make<TH1D>("hist_b_discriminator_ssv_hp" , "B Discriminator for Simple Secondary Vertex High Purity", 100, 0.0, 100.0);
|
| 435 |
hists["hist_ptrel" ] = theDir.make<TH1D>("hist_ptrel" , "P_{T}rel Distribution of Muons", 400, -10, 40);
|
| 436 |
hists["hist_deltaR" ] = theDir.make<TH1D>("hist_deltaR" , "Delta R < 0.7", 20, 0.0, 3.0);
|
| 437 |
hists["hist_eta" ] = theDir.make<TH1D>("hist_eta" , "Eta", 100, -10, 10);
|
| 438 |
//hists["hist_neutralHadronEnergyFraction"] = theDir.make<TH1D>("hist_neutralHadronEnergyFraction", "Neutral Hadron Energy Fraction", 100, -50, 100);
|
| 439 |
// hists["neutralEmEnergyFraction" ] = theDir.make<TH1D>("hist_neutralEmEnergyFraction" , "Neutral Em Energy Fraction", 100, -50, 100);
|
| 440 |
|
| 441 |
hists["hist_nMuons" ] = theDir.make<TH1D>("hist_nMuons" , "Muon Multiplicity", 5, -0.5, 4.5 );
|
| 442 |
// hists["hist_pt_jet0" ] = theDir.make<TH1D>("hist_pt_jet0" , "Pt of the jet with highest pt after matching", 100, -5, 200);
|
| 443 |
//hists["hist_pt_jet1" ] = theDir.make<TH1D>("hist_pt_jet1" , "Pt of the jet with second highest pt after matching", 100, -5, 200);
|
| 444 |
hists["hist_ptmax" ] = theDir.make<TH1D>("hist_ptmax" , "Highest Pt Muon in a jet", 100, 0, 50);
|
| 445 |
hists["hist_ptmax_sec" ] = theDir.make<TH1D>("hist_ptmax_sec" , "Second highest Pt Muon in a jet", 100, 2, 30);
|
| 446 |
cout << "Making functors" << endl;
|
| 447 |
LJMetStudiesSelector theSelector( muonSelParameters,
|
| 448 |
jetSelParameters,
|
| 449 |
metSelParameters,
|
| 450 |
ljmetStudiesParams );
|
| 451 |
|
| 452 |
|
| 453 |
vector<int> const & runs = plotParameters.getParameter<std::vector<int> >("runs");
|
| 454 |
//bool useMC = plotParameters.getParameter<bool>("useMC");
|
| 455 |
|
| 456 |
cout << "About to begin looping" << endl;
|
| 457 |
|
| 458 |
// int counter = 0;
|
| 459 |
// const int size = allMuons.size();
|
| 460 |
// double pt_max[size];
|
| 461 |
|
| 462 |
|
| 463 |
int nev = 0;
|
| 464 |
//loop through each event
|
| 465 |
for (ev.toBegin(); ! ev.atEnd(); ++ev, ++nev) {
|
| 466 |
|
| 467 |
edm::EventBase const & event = ev;
|
| 468 |
|
| 469 |
int run = event.id().run();
|
| 470 |
if ( runs.size() > 0 && find( runs.begin(), runs.end(), run ) == runs.end() ) continue;
|
| 471 |
|
| 472 |
if ( nev % 10000 == 0 ) cout << "Processing run " << event.id().run() << ", event " << event.id().event() << endl;
|
| 473 |
|
| 474 |
std::strbitset ret = theSelector.getBitTemplate();
|
| 475 |
bool passed = theSelector( event, ret );
|
| 476 |
|
| 477 |
|
| 478 |
//
|
| 479 |
//_____ CALO JET _____________________
|
| 480 |
//
|
| 481 |
if ( ret.test("Jet cuts") ) {
|
| 482 |
if ( ret.test("Min jet multiplicity") ) {
|
| 483 |
pt_jet0 = -1000;
|
| 484 |
px_jet0 = -1000;
|
| 485 |
py_jet0 = -1000;
|
| 486 |
pz_jet0 = -1000;
|
| 487 |
phi_jet0 = -1000;
|
| 488 |
eta_jet0 = -1000;
|
| 489 |
pt_jet1 = -1000;
|
| 490 |
px_jet1 = -1000;
|
| 491 |
py_jet1 = -1000;
|
| 492 |
pz_jet1 = -1000;
|
| 493 |
phi_jet1 = -1000;
|
| 494 |
eta_jet1 = -1000;
|
| 495 |
pt_jet2 = -1000;
|
| 496 |
px_jet2 = -1000;
|
| 497 |
py_jet2 = -1000;
|
| 498 |
pz_jet2 = -1000;
|
| 499 |
phi_jet2 = -1000;
|
| 500 |
eta_jet2 = -1000;
|
| 501 |
pt_jet3 = -1000;
|
| 502 |
px_jet3 = -1000;
|
| 503 |
py_jet3 = -1000;
|
| 504 |
pz_jet3 = -1000;
|
| 505 |
phi_jet3 = -1000;
|
| 506 |
eta_jet3 = -1000;
|
| 507 |
pt_jet4 = -1000;
|
| 508 |
px_jet4 = -1000;
|
| 509 |
py_jet4 = -1000;
|
| 510 |
pz_jet4 = -1000;
|
| 511 |
phi_jet4 = -1000;
|
| 512 |
eta_jet4 = -1000;
|
| 513 |
|
| 514 |
bDiscriminator_tche_jet0 = -1000;
|
| 515 |
bDiscriminator_tchp_jet0 = -1000;
|
| 516 |
bDiscriminator_ssv_jet0 = -1000;
|
| 517 |
bDiscriminator_ssv_hp_jet0 = -1000;
|
| 518 |
bDiscriminator_tche_jet1 = -1000;
|
| 519 |
bDiscriminator_tchp_jet1 = -1000;
|
| 520 |
bDiscriminator_ssv_jet1 = -1000;
|
| 521 |
bDiscriminator_ssv_hp_jet1 = -1000;
|
| 522 |
bDiscriminator_tche_jet2 = -1000;
|
| 523 |
bDiscriminator_tchp_jet2 = -1000;
|
| 524 |
bDiscriminator_ssv_jet2 = -1000;
|
| 525 |
bDiscriminator_ssv_hp_jet2 = -1000;
|
| 526 |
bDiscriminator_tche_jet3 = -1000;
|
| 527 |
bDiscriminator_tchp_jet3 = -1000;
|
| 528 |
bDiscriminator_ssv_jet3 = -1000;
|
| 529 |
bDiscriminator_ssv_hp_jet3 = -1000;
|
| 530 |
bDiscriminator_tche_jet4 = -1000;
|
| 531 |
bDiscriminator_tchp_jet4 = -1000;
|
| 532 |
bDiscriminator_ssv_jet4 = -1000;
|
| 533 |
bDiscriminator_ssv_hp_jet4 = -1000;
|
| 534 |
|
| 535 |
|
| 536 |
vector<pat::Jet> const & allJets = theSelector.allJets();
|
| 537 |
int nJets = allJets.size();
|
| 538 |
|
| 539 |
|
| 540 |
|
| 541 |
pat::Jet const & jet0 = theSelector.jet0();
|
| 542 |
|
| 543 |
|
| 544 |
|
| 545 |
pt_jet0 = jet0.pt(); //jet0 is the highest pt jet
|
| 546 |
phi_jet0 = jet0.phi();
|
| 547 |
eta_jet0 = jet0.eta();
|
| 548 |
px_jet0 = jet0.px();
|
| 549 |
py_jet0 = jet0.py();
|
| 550 |
pz_jet0 = jet0.pz();
|
| 551 |
string tagger_tche_jet0 = "trackCountingHighEffBJetTags";
|
| 552 |
string tagger_tchp_jet0 = "trackCountingHighPurBJetTags";
|
| 553 |
string tagger_ssv_jet0 = "simpleSecondaryVertexBJetTags";
|
| 554 |
string tagger_ssv_hp_jet0 = "simpleSecondaryVertexHighPurBJetTags";
|
| 555 |
bDiscriminator_tche_jet0 = jet0.bDiscriminator(tagger_tche_jet0);
|
| 556 |
bDiscriminator_tchp_jet0 = jet0.bDiscriminator(tagger_tchp_jet0);
|
| 557 |
bDiscriminator_ssv_jet0 = jet0.bDiscriminator(tagger_ssv_jet0);
|
| 558 |
bDiscriminator_ssv_hp_jet0 = jet0.bDiscriminator(tagger_ssv_hp_jet0);
|
| 559 |
hists["hist_jet0px"]->Fill( px_jet0 );
|
| 560 |
hists["hist_jet0py"]->Fill( py_jet0 );
|
| 561 |
hists["hist_jet0Pt"]->Fill(pt_jet0);
|
| 562 |
|
| 563 |
|
| 564 |
|
| 565 |
double et0 = jet0.et();
|
| 566 |
|
| 567 |
|
| 568 |
if (nJets > 2 ){
|
| 569 |
|
| 570 |
//For jet 2:
|
| 571 |
pat::Jet const & jet1 = theSelector.jet1();
|
| 572 |
pt_jet1 = jet1.pt();
|
| 573 |
px_jet1 = jet1.px();
|
| 574 |
py_jet1 = jet1.py();
|
| 575 |
pz_jet1 = jet1.pz();
|
| 576 |
eta_jet1 = jet1.eta();
|
| 577 |
phi_jet1 = jet1.phi();
|
| 578 |
string tagger_tche_jet1 = "trackCountingHighEffBJetTags";
|
| 579 |
string tagger_tchp_jet1 = "trackCountingHighPurBJetTags";
|
| 580 |
string tagger_ssv_jet1 = "simpleSecondaryVertexBJetTags";
|
| 581 |
string tagger_ssv_hp_jet1 = "simpleSecondaryVertexHighPurBJetTags";
|
| 582 |
bDiscriminator_tche_jet1 = jet1.bDiscriminator(tagger_tche_jet1);
|
| 583 |
bDiscriminator_tchp_jet1 = jet1.bDiscriminator(tagger_tchp_jet1);
|
| 584 |
bDiscriminator_ssv_jet1 = jet1.bDiscriminator(tagger_ssv_jet1);
|
| 585 |
bDiscriminator_ssv_hp_jet1 = jet1.bDiscriminator(tagger_ssv_hp_jet1);
|
| 586 |
hists["hist_jet1px"]->Fill( px_jet1 );
|
| 587 |
hists["hist_jet1py"]->Fill( py_jet1 );
|
| 588 |
}
|
| 589 |
|
| 590 |
|
| 591 |
|
| 592 |
if (nJets > 3 ){
|
| 593 |
//For jet 3:
|
| 594 |
pat::Jet const & jet2 = theSelector.jet2();
|
| 595 |
pt_jet2 = jet2.pt();
|
| 596 |
px_jet2 = jet2.px();
|
| 597 |
py_jet2 = jet2.py();
|
| 598 |
pz_jet2 = jet2.pz();
|
| 599 |
eta_jet2 = jet2.eta();
|
| 600 |
phi_jet2 = jet2.phi();
|
| 601 |
string tagger_tche_jet2 = "trackCountingHighEffBJetTags";
|
| 602 |
string tagger_tchp_jet2 = "trackCountingHighPurBJetTags";
|
| 603 |
string tagger_ssv_jet2 = "simpleSecondaryVertexBJetTags";
|
| 604 |
string tagger_ssv_hp_jet2 = "simpleSecondaryVertexHighPurBJetTags";
|
| 605 |
bDiscriminator_tche_jet2 = jet2.bDiscriminator(tagger_tche_jet2);
|
| 606 |
bDiscriminator_tchp_jet2 = jet2.bDiscriminator(tagger_tchp_jet2);
|
| 607 |
bDiscriminator_ssv_jet2 = jet2.bDiscriminator(tagger_ssv_jet2);
|
| 608 |
bDiscriminator_ssv_hp_jet2 = jet2.bDiscriminator(tagger_ssv_hp_jet2);
|
| 609 |
|
| 610 |
|
| 611 |
}
|
| 612 |
|
| 613 |
if (nJets > 4){
|
| 614 |
//for jet 4
|
| 615 |
pat::Jet const & jet3 = theSelector.jet3();
|
| 616 |
pt_jet3 = jet3.pt();
|
| 617 |
px_jet3 = jet3.px();
|
| 618 |
py_jet3 = jet3.py();
|
| 619 |
pz_jet3 = jet3.pz();
|
| 620 |
eta_jet3 = jet3.eta();
|
| 621 |
phi_jet3 = jet3.phi();
|
| 622 |
string tagger_tche_jet3 = "trackCountingHighEffBJetTags";
|
| 623 |
string tagger_tchp_jet3 = "trackCountingHighPurBJetTags";
|
| 624 |
string tagger_ssv_jet3 = "simpleSecondaryVertexBJetTags";
|
| 625 |
string tagger_ssv_hp_jet3 = "simpleSecondaryVertexHighPurBJetTags";
|
| 626 |
bDiscriminator_tche_jet3 = jet3.bDiscriminator(tagger_tche_jet3);
|
| 627 |
bDiscriminator_tchp_jet3 = jet3.bDiscriminator(tagger_tchp_jet3);
|
| 628 |
bDiscriminator_ssv_jet3 = jet3.bDiscriminator(tagger_ssv_jet3);
|
| 629 |
bDiscriminator_ssv_hp_jet3 = jet3.bDiscriminator(tagger_ssv_hp_jet3);
|
| 630 |
}
|
| 631 |
|
| 632 |
if (nJets > 5){
|
| 633 |
//for jet 5:
|
| 634 |
pat::Jet const & jet4 = theSelector.jet4();
|
| 635 |
pt_jet4 = jet4.pt();
|
| 636 |
px_jet4 = jet4.px();
|
| 637 |
py_jet4 = jet4.py();
|
| 638 |
pz_jet4 = jet4.pz();
|
| 639 |
eta_jet4 = jet4.eta();
|
| 640 |
phi_jet4 = jet4.phi();
|
| 641 |
string tagger_tche_jet4 = "trackCountingHighEffBJetTags";
|
| 642 |
string tagger_tchp_jet4 = "trackCountingHighPurBJetTags";
|
| 643 |
string tagger_ssv_jet4 = "simpleSecondaryVertexBJetTags";
|
| 644 |
string tagger_ssv_hp_jet4 = "simpleSecondaryVertexHighPurBJetTags";
|
| 645 |
bDiscriminator_tche_jet4 = jet4.bDiscriminator(tagger_tche_jet4);
|
| 646 |
bDiscriminator_tchp_jet4 = jet4.bDiscriminator(tagger_tchp_jet4);
|
| 647 |
bDiscriminator_ssv_jet4 = jet4.bDiscriminator(tagger_ssv_jet4);
|
| 648 |
bDiscriminator_ssv_hp_jet4 = jet4.bDiscriminator(tagger_ssv_hp_jet4);
|
| 649 |
}
|
| 650 |
|
| 651 |
hists["hist_jet1Pt"]->Fill(pt_jet1);
|
| 652 |
|
| 653 |
|
| 654 |
// LorentzVector physicsP4 = jet0.physicsP4();
|
| 655 |
|
| 656 |
|
| 657 |
if (px_jet0 != 0 && py_jet0 != 0 && px_jet1 != 0 && py_jet1 != 0){
|
| 658 |
|
| 659 |
double cos_angle_jets01 = (px_jet0*px_jet1 + py_jet0*py_jet1)/((sqrt(px_jet0*px_jet0 + py_jet0*py_jet0))*(sqrt(px_jet1*px_jet1 +
|
| 660 |
py_jet1*py_jet1)));
|
| 661 |
//std::cout << angle_jets01 << std::endl;
|
| 662 |
|
| 663 |
|
| 664 |
hists["hist_cos_angle_jets01"]->Fill( cos_angle_jets01 );
|
| 665 |
|
| 666 |
}
|
| 667 |
|
| 668 |
|
| 669 |
|
| 670 |
for (vector<pat::Jet>::const_iterator jet = allJets.begin(); //Looping over jets now
|
| 671 |
jet != allJets.end(); jet++){
|
| 672 |
|
| 673 |
|
| 674 |
//std::cout << "Attempt to locate seg fault" << std::endl;
|
| 675 |
|
| 676 |
string tagger1 = "trackCountingHighEffBJetTags";
|
| 677 |
string tagger2 = "softMuonBJetTags";
|
| 678 |
//string tagger3 = "combinedSecondaryVertexBJetTags";
|
| 679 |
string tagger4 = "trackCountingHighPurBJetTags";
|
| 680 |
//string tagger5 = "jetProbabilityBJetTags";
|
| 681 |
string tagger6 = "simpleSecondaryVertexBJetTags";
|
| 682 |
string tagger7 = "simpleSecondaryVertexHighPurBJetTags";
|
| 683 |
|
| 684 |
|
| 685 |
float bDiscriminator1 = jet->bDiscriminator(tagger1);
|
| 686 |
cout << "B-Discriminator for Track Counting Algorithm =" << bDiscriminator1 << endl;
|
| 687 |
|
| 688 |
hists["hist_b_discriminator_tche"]->Fill(bDiscriminator1);
|
| 689 |
|
| 690 |
float bDiscriminator2 = jet->bDiscriminator(tagger2);
|
| 691 |
cout << "B-Discriminator for Soft Lepton Algorithm =" << bDiscriminator2 << endl;
|
| 692 |
|
| 693 |
//float bDiscriminator3 = jet->bDiscriminator(tagger3);
|
| 694 |
//cout << "B-Discriminator for Combined Secondary Vertex Algorithm =" << bDiscriminator3 << endl;
|
| 695 |
|
| 696 |
//hists["hist_b_discriminator_csv"]->Fill(bDiscriminator3);
|
| 697 |
|
| 698 |
float bDiscriminator4 = jet->bDiscriminator(tagger4);
|
| 699 |
|
| 700 |
// std::cout << "SegFault 1" << std::endl;
|
| 701 |
|
| 702 |
hists["hist_b_discriminator_tchp"]->Fill(bDiscriminator4);
|
| 703 |
|
| 704 |
// float bDiscriminator5 = jet->bDiscriminator(tagger5);
|
| 705 |
|
| 706 |
// hists["hist_b_discriminator_jp"]->Fill(bDiscriminator5);
|
| 707 |
|
| 708 |
//cout << "B discriminator for jet probablility" << bDiscriminator5 <<endl;
|
| 709 |
|
| 710 |
float bDiscriminator6 = jet->bDiscriminator(tagger6);
|
| 711 |
|
| 712 |
hists["hist_b_discriminator_ssv"]->Fill(bDiscriminator6);
|
| 713 |
|
| 714 |
float bDiscriminator7 = jet->bDiscriminator(tagger7);
|
| 715 |
|
| 716 |
hists["hist_b_discriminator_ssv_hp"]->Fill(bDiscriminator7);
|
| 717 |
|
| 718 |
string tag_info_label1 = "softMuon";
|
| 719 |
const reco::SoftLeptonTagInfo * tag_info = jet->tagInfoSoftLepton(tag_info_label1);
|
| 720 |
if (tag_info){
|
| 721 |
if (tag_info->leptons()>0){
|
| 722 |
cout << tag_info->properties(0).sip2d << endl ;
|
| 723 |
}
|
| 724 |
}
|
| 725 |
else{
|
| 726 |
// cout << "Tag Info " << tag_info_label1 << "is not found!!!" << endl;
|
| 727 |
}
|
| 728 |
|
| 729 |
string tag_info_label2 = "secondaryVertex";
|
| 730 |
|
| 731 |
std::cout << "SegFault" << std::endl;
|
| 732 |
|
| 733 |
const reco::SecondaryVertexTagInfo * tag_info2 = jet->tagInfoSecondaryVertex(tag_info_label2);
|
| 734 |
if (tag_info2){
|
| 735 |
|
| 736 |
int _n_vertices = jet->tagInfoSecondaryVertex(tag_info_label2)->nVertices();
|
| 737 |
hists["hist_nvertices"]->Fill(_n_vertices);
|
| 738 |
|
| 739 |
}
|
| 740 |
|
| 741 |
|
| 742 |
double eta = jet->eta();
|
| 743 |
|
| 744 |
std::cout << eta << std::endl;
|
| 745 |
|
| 746 |
hists["hist_eta"]->Fill(eta);
|
| 747 |
|
| 748 |
|
| 749 |
|
| 750 |
//double neutralHadronEnergyFraction = jet->neutralHadronEnergyFraction();
|
| 751 |
|
| 752 |
//std::cout << neutralHadronEnergyFraction << std::endl;
|
| 753 |
|
| 754 |
//hists["hist_neutralHadronEnergyFraction"]->Fill(neutralHadronEnergyFraction);
|
| 755 |
|
| 756 |
|
| 757 |
|
| 758 |
/*
|
| 759 |
|
| 760 |
double neutralEmEnergyFraction = jet->neutralEmEnergyFraction ();
|
| 761 |
|
| 762 |
hists["neutralEmEnergyFraction"]->Fill(neutralEmEnergyFraction);
|
| 763 |
|
| 764 |
double chargedHadronEnergyFraction = jet->chargedHadronEnergyFraction();
|
| 765 |
|
| 766 |
hists["chargedHadronEnergyFraction"]->Fill(chargedHadronEnergyFraction);
|
| 767 |
|
| 768 |
double chargedEmEnergyFraction = jet->chargedEmEnergyFraction();
|
| 769 |
|
| 770 |
hists["chargedEmEnergyFraction"]->Fill(chargedEmEnergyFraction);
|
| 771 |
|
| 772 |
double nConstituents = jet->nConstituents();
|
| 773 |
|
| 774 |
hists["nConstituents"]->Fill(nConstituents);
|
| 775 |
|
| 776 |
double chargedMultiplicity = jet->chargedMultiplicity();
|
| 777 |
|
| 778 |
hists["chargedMultiplicity"]->Fill(chargedMultiplicity);
|
| 779 |
|
| 780 |
*/
|
| 781 |
|
| 782 |
|
| 783 |
// double _vx = jet->vx();
|
| 784 |
|
| 785 |
// double _vy = jet->vy();
|
| 786 |
|
| 787 |
// double _vz = jet->vz();
|
| 788 |
|
| 789 |
// cout << "vx = " << _vx << endl;
|
| 790 |
|
| 791 |
// cout << "vy = " << _vy << endl;
|
| 792 |
|
| 793 |
// cout << "vz = " << _vz << endl;
|
| 794 |
|
| 795 |
/*
|
| 796 |
if ( _n_vertices != 0 ){
|
| 797 |
const Vertex a_vertex = jet-> tagInfoSecondaryVertex(tag_info_label2)-> secondaryVertex ( 0 );
|
| 798 |
|
| 799 |
|
| 800 |
if (_n_vertices == 2){
|
| 801 |
|
| 802 |
double _x1 = a_vertex.x();
|
| 803 |
|
| 804 |
double _y1 = a_vertex.y();
|
| 805 |
|
| 806 |
double _z1 = a_vertex.z();
|
| 807 |
|
| 808 |
cout << "x1 = " << _x1 << endl;
|
| 809 |
|
| 810 |
cout << "y1 = " << _y1 << endl;
|
| 811 |
|
| 812 |
cout << "z1 = " << _z1 << endl;
|
| 813 |
|
| 814 |
|
| 815 |
}
|
| 816 |
}
|
| 817 |
|
| 818 |
*/
|
| 819 |
|
| 820 |
|
| 821 |
}
|
| 822 |
|
| 823 |
|
| 824 |
|
| 825 |
/* FIXME: add MC info following this example
|
| 826 |
if ( useMC && jet.genJet() != 0 ) {
|
| 827 |
hists["hist_jetGenEmE"]->Fill( jet.genJet()->emEnergy() );
|
| 828 |
hists["hist_jetGenHadE"]->Fill( jet.genJet()->hadEnergy() );
|
| 829 |
hists["hist_jetEoverGenE"]->Fill( jet.energy() / jet.genJet()->energy() );
|
| 830 |
hists["hist_jetGenEMF"]->Fill( jet.genJet()->emEnergy() / jet.genJet()->energy() );
|
| 831 |
}
|
| 832 |
*/
|
| 833 |
}
|
| 834 |
}
|
| 835 |
|
| 836 |
|
| 837 |
if (ret.test("Muon Cuts") ){
|
| 838 |
if (ret.test("Min muon multiplicity")){
|
| 839 |
|
| 840 |
|
| 841 |
pt_mu0 = -1000;
|
| 842 |
pt_mu1 = -1000;
|
| 843 |
pt_mu2 = -1000;
|
| 844 |
pt_mu3 = -1000;
|
| 845 |
|
| 846 |
pat::Muon const & muon0 = theSelector.muon0();
|
| 847 |
|
| 848 |
|
| 849 |
vector<pat::Muon> const & allMuons = theSelector.allMuons();
|
| 850 |
|
| 851 |
pt_mu0 = muon0.pt();
|
| 852 |
phi_mu0 = muon0.phi();
|
| 853 |
et_mu0 = muon0.et();
|
| 854 |
//int muon_multiplicity = allMuons.size();
|
| 855 |
px = muon0.px();
|
| 856 |
//cout << "muon px = " << px << endl;
|
| 857 |
py = muon0.py();
|
| 858 |
//cout << "muon py = " << py << endl;
|
| 859 |
pz = muon0.pz();
|
| 860 |
//cout << "muon pz = " << pz << endl;
|
| 861 |
nMuons = allMuons.size();
|
| 862 |
cout << "test0" << endl;
|
| 863 |
//TVector3 v(px, py, pz);
|
| 864 |
|
| 865 |
// cout << "v = " << v(1) << "," << v(2) << "," << v(3) << endl;
|
| 866 |
|
| 867 |
hists["hist_muon0Pt"]->Fill(pt_mu0);
|
| 868 |
hists["hist_muon0Phi"]->Fill(phi_mu0);
|
| 869 |
hists["hist_nMuons"]->Fill(nMuons);
|
| 870 |
|
| 871 |
if (nMuons >= 2){
|
| 872 |
|
| 873 |
pat::Muon const & muon1 = theSelector.muon1();
|
| 874 |
pt_mu1 = muon1.pt();
|
| 875 |
|
| 876 |
}
|
| 877 |
|
| 878 |
if (nMuons >= 3){
|
| 879 |
|
| 880 |
pat::Muon const & muon2 = theSelector.muon2();
|
| 881 |
pt_mu2 = muon2.pt();
|
| 882 |
|
| 883 |
}
|
| 884 |
|
| 885 |
if (nMuons >= 4){
|
| 886 |
|
| 887 |
pat::Muon const & muon3 = theSelector.muon3();
|
| 888 |
pt_mu3 = muon3.pt();
|
| 889 |
|
| 890 |
}
|
| 891 |
|
| 892 |
|
| 893 |
|
| 894 |
|
| 895 |
|
| 896 |
}
|
| 897 |
}
|
| 898 |
|
| 899 |
|
| 900 |
// cout << "px =" << px << endl;
|
| 901 |
|
| 902 |
double px_mu_jet = px_jet0 + px;
|
| 903 |
double py_mu_jet = py_jet0 + py;
|
| 904 |
double pz_mu_jet = pz_jet0 + pz;
|
| 905 |
|
| 906 |
vector<pat::Jet> const & allJets = theSelector.allJets();
|
| 907 |
|
| 908 |
vector<pat::Muon> const & allMuons = theSelector.allMuons();
|
| 909 |
|
| 910 |
//pat::Jet const & jet0 = theSelector.jet0();
|
| 911 |
//pat::Jet const & jet1 = theSelector.jet1();
|
| 912 |
|
| 913 |
int mu = 0;
|
| 914 |
|
| 915 |
//int counterMu = 0;
|
| 916 |
//const int sizeMu = allMuons.size();
|
| 917 |
//double pt_max[sizeMu];
|
| 918 |
vector<double> pt_max;
|
| 919 |
vector<double> pt_max_sec;
|
| 920 |
|
| 921 |
|
| 922 |
|
| 923 |
|
| 924 |
for (vector<pat::Jet>::const_iterator jet = allJets.begin(); //Looping over jets now
|
| 925 |
jet != allJets.end(); jet++){
|
| 926 |
|
| 927 |
//int counter = 0;
|
| 928 |
//const int sizeMu = allMuons.size();
|
| 929 |
//double pt_max[sizeMu];
|
| 930 |
|
| 931 |
for ( vector<pat::Muon>::const_iterator mu = allMuons.begin();
|
| 932 |
mu != allMuons.end(); mu++){
|
| 933 |
|
| 934 |
//for(i=0; i!= allMuons.size()-1; i++){
|
| 935 |
double px_allmu = mu->px();
|
| 936 |
double py_allmu = mu->py();
|
| 937 |
double pz_allmu = mu->pz();
|
| 938 |
|
| 939 |
double px_alljets = jet->px();
|
| 940 |
double py_alljets = jet->py();
|
| 941 |
double pz_alljets = jet->pz();
|
| 942 |
|
| 943 |
double px_amu_ajet = px_allmu + px_alljets;
|
| 944 |
double py_amu_ajet = py_allmu + py_alljets;
|
| 945 |
double pz_amu_ajet = pz_allmu + pz_alljets;
|
| 946 |
|
| 947 |
|
| 948 |
TVector3 v4(px_allmu, py_allmu, pz_allmu);
|
| 949 |
|
| 950 |
TVector3 v5(px_alljets, py_alljets, pz_alljets);
|
| 951 |
|
| 952 |
double deltaR_mu_jet = v4.DeltaR(v5);
|
| 953 |
|
| 954 |
if (deltaR_mu_jet < 0.5){
|
| 955 |
|
| 956 |
if (px_amu_ajet != 0 && py_amu_ajet != 0 && pz_amu_ajet != 0){
|
| 957 |
|
| 958 |
|
| 959 |
double ptrel = sqrt((py_allmu*pz_amu_ajet - pz_allmu*py_amu_ajet)*(py_allmu*pz_amu_ajet - pz_allmu*py_amu_ajet) + (pz_allmu*px_amu_ajet -
|
| 960 |
px_allmu*pz_amu_ajet)*(pz_allmu*px_amu_ajet - px_allmu*pz_amu_ajet) + (px_allmu*pz_amu_ajet - py_allmu*px_amu_ajet)*(px_allmu*pz_amu_ajet - py_allmu*px_amu_ajet))/(sqrt(px_amu_ajet*px_amu_ajet + py_amu_ajet*py_amu_ajet + pz_amu_ajet*pz_amu_ajet));
|
| 961 |
|
| 962 |
// std::cout << ptrel << std::endl;
|
| 963 |
|
| 964 |
hists["hist_ptrel"]->Fill(ptrel); }
|
| 965 |
|
| 966 |
|
| 967 |
pt_max.push_back(mu->pt());
|
| 968 |
cout << pt_max[0] << endl;
|
| 969 |
hists["hist_ptmax"]->Fill(pt_max[0]);
|
| 970 |
if (allMuons.size() >= 2 && pt_max[1] != 0.0){
|
| 971 |
|
| 972 |
hists["hist_ptmax_sec"]->Fill(pt_max[1]);
|
| 973 |
}
|
| 974 |
|
| 975 |
|
| 976 |
}
|
| 977 |
|
| 978 |
}
|
| 979 |
|
| 980 |
}
|
| 981 |
|
| 982 |
if (ret.test("Min muon multiplicity") && px !=0 && py !=0 && pz !=0 && px_jet1 != 0 && py_jet1 != 0 && pz_jet1 !=0 ){
|
| 983 |
|
| 984 |
TVector3 v1(px_jet1, py_jet1, pz_jet1);
|
| 985 |
|
| 986 |
// cout << "TVector3 v1 = " << v1(1) << "," << v1(2) << "," << v1(3) << endl;
|
| 987 |
|
| 988 |
TVector3 v2(px, py, pz);
|
| 989 |
|
| 990 |
cout << "TVector3 v2 = " << v2(0) << "," << v2(1) << "," << v2(2) << endl;
|
| 991 |
|
| 992 |
|
| 993 |
double deltaR = v1.DeltaR(v2);
|
| 994 |
|
| 995 |
if (deltaR < 0.5 ){
|
| 996 |
|
| 997 |
cout << "deltaR is less than 0.5" << endl;
|
| 998 |
|
| 999 |
hists["hist_deltaR"]->Fill(deltaR);
|
| 1000 |
|
| 1001 |
}
|
| 1002 |
|
| 1003 |
|
| 1004 |
}
|
| 1005 |
|
| 1006 |
tree["Kinematic_Variables"]->Fill();
|
| 1007 |
|
| 1008 |
|
| 1009 |
} // end loop over events
|
| 1010 |
|
| 1011 |
cout << "Selection" << endl;
|
| 1012 |
theSelector.print(std::cout);
|
| 1013 |
|
| 1014 |
|
| 1015 |
|
| 1016 |
|
| 1017 |
return 0;
|
| 1018 |
}
|
| 1019 |
|
