27 #include "TDirectory.h"
32 #include <Math/VectorUtil.h>
81 if (objectType_ ==
"electron")
83 else if (objectType_ ==
"muon")
85 else if (objectType_ ==
"lJets" || objectType_ ==
"bJets")
87 else if (objectType_ ==
"met")
89 else if (objectType_ ==
"tau")
91 if (objectType_ !=
"met") {
109 std::vector<reco::Particle *> p4gen, p4rec;
114 if (genEvt->particles().size() < 10)
121 for (
size_t e = 0;
e < electrons->size();
e++) {
122 for (
size_t p = 0;
p < genEvt->particles().size();
p++) {
123 if ((
std::abs(genEvt->particles()[
p].pdgId()) == 11) &&
133 for (
size_t m = 0;
m < muons->size();
m++) {
134 for (
size_t p = 0;
p < genEvt->particles().size();
p++) {
135 if ((
std::abs(genEvt->particles()[
p].pdgId()) == 13) &&
145 if (jets->size() >= 4) {
146 for (
unsigned int j = 0;
j < 4;
j++) {
147 for (
size_t p = 0;
p < genEvt->particles().size();
p++) {
148 if ((
std::abs(genEvt->particles()[
p].pdgId()) < 5) &&
159 if (jets->size() >= 4) {
160 for (
unsigned int j = 0;
j < 4;
j++) {
161 for (
size_t p = 0;
p < genEvt->particles().size();
p++) {
162 if ((
std::abs(genEvt->particles()[
p].pdgId()) == 5) &&
173 if (!mets->empty()) {
174 if (genEvt->isSemiLeptonic() && genEvt->singleNeutrino() !=
nullptr &&
183 for (std::vector<pat::Tau>::const_iterator
tau = taus->begin();
tau != taus->end(); ++
tau) {
196 for (
unsigned m = 0;
m < p4gen.size();
m++) {
197 double Egen = p4gen[
m]->
energy();
198 double Thetagen = p4gen[
m]->theta();
199 double Phigen = p4gen[
m]->phi();
200 double Etgen = p4gen[
m]->et();
201 double Etagen = p4gen[
m]->eta();
202 double Ecal = p4rec[
m]->energy();
203 double Thetacal = p4rec[
m]->theta();
204 double Phical = p4rec[
m]->phi();
205 double Etcal = p4rec[
m]->et();
206 double Etacal = p4rec[
m]->eta();
207 double phidiff = Phical - Phigen;
208 if (phidiff > 3.14159)
209 phidiff = 2. * 3.14159 -
phidiff;
210 if (phidiff < -3.14159)
211 phidiff = -phidiff - 2. * 3.14159;
239 ROOT::Math::SVector<double, 3> pcalvec(p4rec[
m]->px(), p4rec[
m]->py(), p4rec[
m]->pz());
240 ROOT::Math::SVector<double, 3> pgenvec(p4gen[
m]->px(), p4gen[
m]->py(), p4gen[
m]->pz());
242 ROOT::Math::SVector<double, 3> u_z(0, 0, 1);
243 ROOT::Math::SVector<double, 3> u_1 = ROOT::Math::Unit(pcalvec);
244 ROOT::Math::SVector<double, 3> u_3 = ROOT::Math::Cross(u_z, u_1) / ROOT::Math::Mag(ROOT::Math::Cross(u_z, u_1));
245 ROOT::Math::SVector<double, 3> u_2 = ROOT::Math::Cross(u_1, u_3) / ROOT::Math::Mag(ROOT::Math::Cross(u_1, u_3));
250 double agen = ROOT::Math::Dot(pgenvec, u_1) / ROOT::Math::Mag(pcalvec);
251 double bgen = ROOT::Math::Dot(pgenvec, u_2);
252 double cgen = ROOT::Math::Dot(pgenvec, u_3);
253 double dgen = Egen / Ecal;
261 hResPtEtaBin[4][etabin][ptbin]->Fill(Thetacal - Thetagen);
278 TString resObsName[8] = {
"_ares",
"_bres",
"_cres",
"_dres",
"_thres",
"_phres",
"_etres",
"_etares"};
279 int resObsNrBins = 120;
282 std::vector<double> resObsMin, resObsMax;
284 resObsMin.push_back(-0.15);
285 resObsMin.push_back(-0.2);
286 resObsMin.push_back(-0.1);
287 resObsMin.push_back(-0.15);
288 resObsMin.push_back(-0.0012);
289 resObsMin.push_back(-0.009);
290 resObsMin.push_back(-16);
291 resObsMin.push_back(-0.0012);
292 resObsMax.push_back(0.15);
293 resObsMax.push_back(0.2);
294 resObsMax.push_back(0.1);
295 resObsMax.push_back(0.15);
296 resObsMax.push_back(0.0012);
297 resObsMax.push_back(0.009);
298 resObsMax.push_back(16);
299 resObsMax.push_back(0.0012);
301 resObsMin.push_back(-0.15);
302 resObsMin.push_back(-0.1);
303 resObsMin.push_back(-0.05);
304 resObsMin.push_back(-0.15);
305 resObsMin.push_back(-0.004);
306 resObsMin.push_back(-0.003);
307 resObsMin.push_back(-8);
308 resObsMin.push_back(-0.004);
309 resObsMax.push_back(0.15);
310 resObsMax.push_back(0.1);
311 resObsMax.push_back(0.05);
312 resObsMax.push_back(0.15);
313 resObsMax.push_back(0.004);
314 resObsMax.push_back(0.003);
315 resObsMax.push_back(8);
316 resObsMax.push_back(0.004);
318 resObsMin.push_back(-1.);
319 resObsMin.push_back(-10.);
320 resObsMin.push_back(-10);
321 resObsMin.push_back(-1.);
322 resObsMin.push_back(-0.1);
323 resObsMin.push_back(-0.1);
324 resObsMin.push_back(-80);
325 resObsMin.push_back(-0.1);
326 resObsMax.push_back(1.);
327 resObsMax.push_back(10.);
328 resObsMax.push_back(10);
329 resObsMax.push_back(1.);
330 resObsMax.push_back(0.1);
331 resObsMax.push_back(0.1);
332 resObsMax.push_back(50);
333 resObsMax.push_back(0.1);
335 resObsMin.push_back(-1.);
336 resObsMin.push_back(-10.);
337 resObsMin.push_back(-10.);
338 resObsMin.push_back(-1.);
339 resObsMin.push_back(-0.4);
340 resObsMin.push_back(-0.6);
341 resObsMin.push_back(-80);
342 resObsMin.push_back(-0.6);
343 resObsMax.push_back(1.);
344 resObsMax.push_back(10.);
345 resObsMax.push_back(10.);
346 resObsMax.push_back(1.);
347 resObsMax.push_back(0.4);
348 resObsMax.push_back(0.6);
349 resObsMax.push_back(80);
350 resObsMax.push_back(0.6);
352 resObsMin.push_back(-2.);
353 resObsMin.push_back(-150.);
354 resObsMin.push_back(-150.);
355 resObsMin.push_back(-2.);
356 resObsMin.push_back(-6);
357 resObsMin.push_back(-6);
358 resObsMin.push_back(-180);
359 resObsMin.push_back(-6);
360 resObsMax.push_back(3.);
361 resObsMax.push_back(150.);
362 resObsMax.push_back(150.);
363 resObsMax.push_back(3.);
364 resObsMax.push_back(6);
365 resObsMax.push_back(6);
366 resObsMax.push_back(180);
367 resObsMax.push_back(6);
370 const char *resObsVsPtFit[8] = {
"[0]+[1]*exp(-[2]*x)",
371 "[0]+[1]*exp(-[2]*x)",
372 "[0]+[1]*exp(-[2]*x)",
373 "[0]+[1]*exp(-[2]*x)",
374 "[0]+[1]*exp(-[2]*x)",
375 "[0]+[1]*exp(-[2]*x)",
377 "[0]+[1]*exp(-[2]*x)"};
391 etabins =
new double[2];
397 for (Int_t ro = 0; ro < 8; ro++) {
398 for (Int_t etab = 0; etab <
etanrbins; etab++) {
399 for (Int_t ptb = 0; ptb <
ptnrbins; ptb++) {
401 obsName += resObsName[ro];
402 obsName +=
"_etabin";
406 hResPtEtaBin[ro][etab][ptb] = fs->
make<TH1F>(obsName, obsName, resObsNrBins, resObsMin[ro], resObsMax[ro]);
410 obsName2 += resObsName[ro];
411 obsName2 +=
"_etabin";
415 fs->
make<TF1>(
"F_" + obsName2, resObsVsPtFit[ro],
pTbinVals_[0], pTbinVals_[pTbinVals_.size() - 1]);
418 tResVar = fs->
make<TTree>(
"tResVar",
"Resolution tree");
426 TString resObsName2[8] = {
"a",
"b",
"c",
"d",
"theta",
"phi",
"et",
"eta"};
432 tResVar->Branch(
"Pt", &pt,
"Pt/D");
433 tResVar->Branch(
"Eta", &eta,
"Eta/D");
434 tResVar->Branch(
"ro", &ro,
"ro/I");
435 tResVar->Branch(
"value", &value,
"value/D");
436 tResVar->Branch(
"error", &error,
"error/D");
438 for (ro = 0; ro < 8; ro++) {
439 for (
int etab = 0; etab <
etanrbins; etab++) {
442 for (
int ptb = 0; ptb <
ptnrbins; ptb++) {
445 double maxcontent = 0.;
447 for (
int nb = 1; nb <
hResPtEtaBin[ro][etab][ptb]->GetNbinsX(); nb++) {
448 if (
hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb) > maxcontent) {
449 maxcontent =
hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb);
455 hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin + range));
481 edm::LogProblem(
"SummaryError") <<
"No plots filled for light jets \n";
494 edm::LogVerbatim(
"MainResults") <<
" ----------------------------------------------";
495 edm::LogVerbatim(
"MainResults") <<
" ----------------------------------------------";
497 edm::LogVerbatim(
"MainResults") <<
" ----------------------------------------------";
498 edm::LogVerbatim(
"MainResults") <<
" ----------------------------------------------";
500 for (ro = 0; ro < 8; ro++) {
502 edm::LogVerbatim(
"MainResults") <<
"\n Resolutions on " << resObsName2[ro] <<
"\n";
504 for (
int etab = 0; etab <
etanrbins; etab++) {
505 if (nrFilled != 0 && ro != 6) {
509 <<
"+" <<
fResEtaBin[ro][etab]->GetParameter(1) <<
"*exp(-(" <<
fResEtaBin[ro][etab]->GetParameter(2)
513 <<
"else if(fabs(eta)<" <<
etabinVals_[etab + 1] <<
") res = " <<
fResEtaBin[ro][etab]->GetParameter(0)
514 <<
"+" <<
fResEtaBin[ro][etab]->GetParameter(1) <<
"*exp(-(" <<
fResEtaBin[ro][etab]->GetParameter(2)
517 }
else if (nrFilled != 0 && ro == 6) {
521 <<
"+" <<
fResEtaBin[ro][etab]->GetParameter(1) <<
"*pt;";
524 <<
"else if(fabs(eta)<" <<
etabinVals_[etab + 1] <<
") res = " <<
fResEtaBin[ro][etab]->GetParameter(0)
525 <<
"+" <<
fResEtaBin[ro][etab]->GetParameter(1) <<
"*pt;";
530 }
else if (nrFilled != 0 &&
objectType_ ==
"met") {
531 for (ro = 0; ro < 8; ro++) {
533 edm::LogVerbatim(
"MainResults") <<
"\n Resolutions on " << resObsName2[ro] <<
"\n";
535 if (nrFilled != 0 && ro != 6) {
537 <<
fResEtaBin[ro][0]->GetParameter(1) <<
"*exp(-("
538 <<
fResEtaBin[ro][0]->GetParameter(2) <<
"*pt));";
539 }
else if (nrFilled != 0 && ro == 6) {
541 <<
fResEtaBin[ro][0]->GetParameter(1) <<
"*pt;";
Log< level::Info, true > LogVerbatim
const Candidate * mother(size_type=0) const override
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
int status() const final
status word
T * make(const Args &...args) const
make new ROOT object
double phidiff(double phi)
Normalized difference in azimuthal angles to a range between .
edm::EDGetTokenT< TtGenEvent > genEvtToken_
size_t numberOfMothers() const override
number of mothers
const LorentzVector & p4() const final
four-momentum Lorentz vector
const uint16_t range(const Frame &aFrame)
int pdgId() const final
PDG identifier.
edm::EDGetTokenT< std::vector< pat::Jet > > jetsToken_
edm::EDGetTokenT< std::vector< pat::Tau > > tausToken_
std::vector< double > etabinVals_
~ResolutionCreator() override
Abs< T >::type abs(const T &t)
TH1F * hResEtaBin[10][20]
TF1 * fResPtEtaBin[10][20][20]
virtual int pdgId() const =0
PDG identifier.
TH1F * hResPtEtaBin[10][20][20]
void analyze(const edm::Event &, const edm::EventSetup &) override
ResolutionCreator(const edm::ParameterSet &)
edm::EDGetTokenT< std::vector< pat::Electron > > electronsToken_
T getParameter(std::string const &) const
edm::EDGetTokenT< std::vector< pat::MET > > metsToken_
std::vector< double > pTbinVals_
edm::EDGetTokenT< std::vector< pat::Muon > > muonsToken_
Log< level::Error, true > LogProblem
double energy() const final
energy