27 #include "TDirectory.h"
32 #include <Math/VectorUtil.h>
102 std::vector<reco::Particle *> p4gen, p4rec;
107 if(genEvt->particles().size()<10)
return;
112 for(
size_t e=0;
e<electrons->size();
e++) {
113 for(
size_t p=0;
p<genEvt->particles().size();
p++){
114 if( (
std::abs(genEvt->particles()[
p].pdgId()) == 11) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[
p].p4(), (*electrons)[
e].p4()) <
minDR_) ) {
124 for(
size_t m=0;
m<muons->size();
m++) {
125 for(
size_t p=0;
p<genEvt->particles().size();
p++){
126 if( (
std::abs(genEvt->particles()[
p].pdgId()) == 13) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[
p].p4(), (*muons)[
m].p4()) <
minDR_) ) {
136 if(jets->size()>=4) {
137 for(
unsigned int j = 0;
j<4;
j++){
138 for(
size_t p=0;
p<genEvt->particles().size();
p++){
139 if( (
std::abs(genEvt->particles()[
p].pdgId()) < 5) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[
p].p4(), (*jets)[
j].p4())<
minDR_) ){
150 if(jets->size()>=4) {
151 for(
unsigned int j = 0;
j<4;
j++){
152 for(
size_t p=0;
p<genEvt->particles().size();
p++){
153 if( (
std::abs(genEvt->particles()[
p].pdgId()) == 5) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[
p].p4(), (*jets)[
j].p4())<
minDR_) ) {
164 if(mets->size()>=1) {
165 if( genEvt->isSemiLeptonic() && genEvt->singleNeutrino() != 0 && ROOT::Math::VectorUtil::DeltaR(genEvt->singleNeutrino()->p4(), (*mets)[0].p4()) <
minDR_) {
174 for(std::vector<pat::Tau>::const_iterator
tau = taus->begin();
tau != taus->end(); ++
tau) {
182 ROOT::Math::VectorUtil::DeltaR(genLepton.
p4(),
tau->p4()) <
minDR_ ) {
189 for(
unsigned m=0;
m<p4gen.size();
m++){
190 double Egen = p4gen[
m]->
energy();
191 double Thetagen = p4gen[
m]->theta();
192 double Phigen = p4gen[
m]->phi();
193 double Etgen = p4gen[
m]->et();
194 double Etagen = p4gen[
m]->eta();
195 double Ecal = p4rec[
m]->energy();
196 double Thetacal = p4rec[
m]->theta();
197 double Phical = p4rec[
m]->phi();
198 double Etcal = p4rec[
m]->et();
199 double Etacal = p4rec[
m]->eta();
200 double phidiff = Phical- Phigen;
201 if(phidiff>3.14159) phidiff = 2.*3.14159 -
phidiff;
202 if(phidiff<-3.14159) phidiff = -phidiff - 2.*3.14159;
228 ROOT::Math::SVector<double,3> pcalvec(p4rec[
m]->px(),p4rec[
m]->py(),p4rec[
m]->pz());
229 ROOT::Math::SVector<double,3> pgenvec(p4gen[
m]->px(),p4gen[
m]->py(),p4gen[
m]->pz());
231 ROOT::Math::SVector<double,3> u_z(0,0,1);
232 ROOT::Math::SVector<double,3> u_1 = ROOT::Math::Unit(pcalvec);
233 ROOT::Math::SVector<double,3> u_3 = ROOT::Math::Cross(u_z,u_1)/ROOT::Math::Mag(ROOT::Math::Cross(u_z,u_1));
234 ROOT::Math::SVector<double,3> u_2 = ROOT::Math::Cross(u_1,u_3)/ROOT::Math::Mag(ROOT::Math::Cross(u_1,u_3));
239 double agen = ROOT::Math::Dot(pgenvec,u_1)/ROOT::Math::Mag(pcalvec);
240 double bgen = ROOT::Math::Dot(pgenvec,u_2);
241 double cgen = ROOT::Math::Dot(pgenvec,u_3);
242 double dgen = Egen/Ecal;
270 TString resObsName[8] = {
"_ares",
"_bres",
"_cres",
"_dres",
"_thres",
"_phres",
"_etres",
"_etares"};
271 int resObsNrBins = 120;
273 std::vector<double> resObsMin, resObsMax;
275 resObsMin.push_back(-0.15); resObsMin.push_back(-0.2); resObsMin.push_back(-0.1); resObsMin.push_back(-0.15); resObsMin.push_back(-0.0012); resObsMin.push_back(-0.009); resObsMin.push_back(-16); resObsMin.push_back(-0.0012);
276 resObsMax.push_back( 0.15); resObsMax.push_back( 0.2); resObsMax.push_back( 0.1); resObsMax.push_back( 0.15); resObsMax.push_back( 0.0012); resObsMax.push_back( 0.009); resObsMax.push_back( 16); resObsMax.push_back( 0.0012);
278 resObsMin.push_back(-0.15); resObsMin.push_back(-0.1); resObsMin.push_back(-0.05); resObsMin.push_back(-0.15); resObsMin.push_back(-0.004); resObsMin.push_back(-0.003); resObsMin.push_back(-8); resObsMin.push_back(-0.004);
279 resObsMax.push_back( 0.15); resObsMax.push_back( 0.1); resObsMax.push_back( 0.05); resObsMax.push_back( 0.15); resObsMax.push_back( 0.004); resObsMax.push_back( 0.003); resObsMax.push_back( 8); resObsMax.push_back( 0.004);
281 resObsMin.push_back(-1.); resObsMin.push_back(-10.); resObsMin.push_back(-10); resObsMin.push_back(-1.); resObsMin.push_back(-0.1); resObsMin.push_back(-0.1); resObsMin.push_back(-80); resObsMin.push_back(-0.1);
282 resObsMax.push_back( 1.); resObsMax.push_back( 10.); resObsMax.push_back( 10); resObsMax.push_back( 1.); resObsMax.push_back( 0.1); resObsMax.push_back( 0.1); resObsMax.push_back( 50); resObsMax.push_back( 0.1);
284 resObsMin.push_back(-1.); resObsMin.push_back(-10.); resObsMin.push_back(-10.); resObsMin.push_back(-1.); resObsMin.push_back(-0.4); resObsMin.push_back(-0.6); resObsMin.push_back( -80); resObsMin.push_back(-0.6);
285 resObsMax.push_back( 1.); resObsMax.push_back( 10.); resObsMax.push_back( 10.); resObsMax.push_back( 1.); resObsMax.push_back( 0.4); resObsMax.push_back( 0.6); resObsMax.push_back( 80); resObsMax.push_back( 0.6);
287 resObsMin.push_back(-2.); resObsMin.push_back(-150.); resObsMin.push_back(-150.); resObsMin.push_back(-2.); resObsMin.push_back(-6); resObsMin.push_back(-6); resObsMin.push_back( -180); resObsMin.push_back(-6);
288 resObsMax.push_back( 3.); resObsMax.push_back( 150.); resObsMax.push_back( 150.); resObsMax.push_back( 3.); resObsMax.push_back( 6); resObsMax.push_back( 6); resObsMax.push_back( 180); resObsMax.push_back( 6);
291 const char* resObsVsPtFit[8] = {
"[0]+[1]*exp(-[2]*x)",
292 "[0]+[1]*exp(-[2]*x)",
293 "[0]+[1]*exp(-[2]*x)",
294 "[0]+[1]*exp(-[2]*x)",
295 "[0]+[1]*exp(-[2]*x)",
296 "[0]+[1]*exp(-[2]*x)",
298 "[0]+[1]*exp(-[2]*x)"
302 double *ptbins =
new double[
pTbinVals_.size()];
311 etabins =
new double[2];
312 etabins[0] = 0; etabins[1] = 5.;
317 for(Int_t ro=0; ro<8; ro++) {
318 for(Int_t etab=0; etab<
etanrbins; etab++) {
319 for(Int_t ptb=0; ptb<
ptnrbins; ptb++) {
320 TString obsName =
objectType_; obsName += resObsName[ro]; obsName +=
"_etabin"; obsName += etab; obsName +=
"_ptbin";
322 hResPtEtaBin[ro][etab][ptb] = fs->
make<TH1F>(obsName,obsName,resObsNrBins,resObsMin[ro],resObsMax[ro]);
325 TString obsName2 =
objectType_; obsName2 += resObsName[ro]; obsName2 +=
"_etabin"; obsName2 += etab;
327 fResEtaBin[ro][etab] = fs->
make<TF1>(
"F_"+obsName2,resObsVsPtFit[ro],
pTbinVals_[0],pTbinVals_[pTbinVals_.size()-1]);
330 tResVar = fs->
make< TTree >(
"tResVar",
"Resolution tree");
340 TString resObsName2[8] = {
"a",
"b",
"c",
"d",
"theta",
"phi",
"et",
"eta"};
346 tResVar->Branch(
"Pt",&pt,
"Pt/D");
347 tResVar->Branch(
"Eta",&eta,
"Eta/D");
348 tResVar->Branch(
"ro",&ro,
"ro/I");
349 tResVar->Branch(
"value",&value,
"value/D");
350 tResVar->Branch(
"error",&error,
"error/D");
352 for(ro=0; ro<8; ro++) {
353 for(
int etab=0; etab<
etanrbins; etab++) {
356 for(
int ptb=0; ptb<
ptnrbins; ptb++) {
359 double maxcontent = 0.;
361 for(
int nb=1; nb<
hResPtEtaBin[ro][etab][ptb]->GetNbinsX(); nb ++){
362 if (
hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb)>maxcontent) {
363 maxcontent =
hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb);
367 int range = (int)(
hResPtEtaBin[ro][etab][ptb]->GetNbinsX()/6);
369 hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin+range));
399 edm::LogVerbatim (
"MainResults") <<
" ----------------------------------------------";
400 edm::LogVerbatim (
"MainResults") <<
" ----------------------------------------------";
402 edm::LogVerbatim (
"MainResults") <<
" ----------------------------------------------";
403 edm::LogVerbatim (
"MainResults") <<
" ----------------------------------------------";
405 for(ro=0; ro<8; ro++) {
407 edm::LogVerbatim (
"MainResults") <<
"\n Resolutions on " << resObsName2[ro] <<
"\n";
409 for(
int etab=0; etab<
etanrbins; etab++) {
410 if(nrFilled != 0 && ro != 6) {
414 <<
"*exp(-(" <<
fResEtaBin[ro][etab]->GetParameter(2) <<
"*pt));";
418 <<
"*exp(-(" <<
fResEtaBin[ro][etab]->GetParameter(2) <<
"*pt));";
420 }
else if(nrFilled != 0 && ro == 6){
435 for(ro=0; ro<8; ro++) {
437 edm::LogVerbatim (
"MainResults") <<
"\n Resolutions on " << resObsName2[ro] <<
"\n";
439 if(nrFilled != 0 && ro != 6) {
442 <<
"*exp(-(" <<
fResEtaBin[ro][0]->GetParameter(2) <<
"*pt));";
443 }
else if(nrFilled != 0 && ro == 6){
virtual double energy() const GCC11_FINAL
energy
T getParameter(std::string const &) const
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
virtual const LorentzVector & p4() const GCC11_FINAL
four-momentum Lorentz vector
virtual int pdgId() const GCC11_FINAL
PDG identifier.
#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)
double phidiff(double phi)
Normalized difference in azimuthal angles to a range between .
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::vector< double > etabinVals_
virtual int status() const GCC11_FINAL
status word
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
virtual size_t numberOfMothers() const
number of mothers
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
TH1F * hResEtaBin[10][20]
TF1 * fResPtEtaBin[10][20][20]
virtual int pdgId() const =0
PDG identifier.
TH1F * hResPtEtaBin[10][20][20]
ResolutionCreator(const edm::ParameterSet &)
T * make() const
make new ROOT object
std::vector< double > pTbinVals_
virtual const Candidate * mother(size_type=0) const
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...