27 #include "TDirectory.h"
32 #include <Math/VectorUtil.h>
109 std::vector<reco::Particle *> p4gen, p4rec;
114 if (
genEvt->particles().size() < 10)
122 for (
size_t p = 0;
p <
genEvt->particles().size();
p++) {
133 for (
size_t m = 0;
m <
muons->size();
m++) {
134 for (
size_t p = 0;
p <
genEvt->particles().size();
p++) {
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++) {
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++) {
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;
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)"};
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";
418 tResVar = fs->
make<TTree>(
"tResVar",
"Resolution tree");
426 TString resObsName2[8] = {
"a",
"b",
"c",
"d",
"theta",
"phi",
"et",
"eta"};
434 tResVar->Branch(
"ro", &ro,
"ro/I");
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);
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++) {
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;";
531 for (ro = 0; ro < 8; ro++) {
533 edm::LogVerbatim(
"MainResults") <<
"\n Resolutions on " << resObsName2[ro] <<
"\n";
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;";