Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:24 

0001 // system include files
0002 #include <memory>
0003 #include <string>
0004 
0005 // user include files
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 //needed for TFileService
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 //needed for MessageLogger
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
0020 
0021 #include "DataFormats/PatCandidates/interface/Muon.h"
0022 #include "DataFormats/PatCandidates/interface/Jet.h"
0023 #include "DataFormats/PatCandidates/interface/MET.h"
0024 #include "DataFormats/PatCandidates/interface/Tau.h"
0025 #include "DataFormats/PatCandidates/interface/Electron.h"
0026 
0027 #include "TDirectory.h"
0028 #include "TH1F.h"
0029 #include "TF1.h"
0030 #include "TTree.h"
0031 
0032 #include <Math/VectorUtil.h>
0033 
0034 //
0035 // class declaration
0036 //
0037 
0038 class ResolutionCreator : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0039 public:
0040   explicit ResolutionCreator(const edm::ParameterSet &);
0041 
0042 private:
0043   void beginJob() override;
0044   void analyze(const edm::Event &, const edm::EventSetup &) override;
0045   void endJob() override;
0046 
0047   // ----------member data ---------------------------
0048   edm::EDGetTokenT<TtGenEvent> genEvtToken_;
0049   std::string objectType_, labelName_;
0050   edm::EDGetTokenT<std::vector<pat::Electron> > electronsToken_;
0051   edm::EDGetTokenT<std::vector<pat::Muon> > muonsToken_;
0052   edm::EDGetTokenT<std::vector<pat::Jet> > jetsToken_;
0053   edm::EDGetTokenT<std::vector<pat::MET> > metsToken_;
0054   edm::EDGetTokenT<std::vector<pat::Tau> > tausToken_;
0055   std::vector<double> etabinVals_, pTbinVals_;
0056   double minDR_;
0057   int matchingAlgo_;
0058   bool useMaxDist_;
0059   bool useDeltaR_;
0060   double maxDist_;
0061   int ptnrbins, etanrbins;
0062   int nrFilled;
0063 
0064   //Histograms are booked in the beginJob() method
0065   TF1 *fResPtEtaBin[10][20][20];
0066   TF1 *fResEtaBin[10][20];
0067   TH1F *hResPtEtaBin[10][20][20];
0068   TH1F *hResEtaBin[10][20];
0069   TTree *tResVar;
0070 };
0071 
0072 //
0073 // constructors and destructor
0074 //
0075 ResolutionCreator::ResolutionCreator(const edm::ParameterSet &iConfig) {
0076   usesResource("TFileService");
0077 
0078   // input parameters
0079   genEvtToken_ = consumes<TtGenEvent>(edm::InputTag("genEvt"));
0080   objectType_ = iConfig.getParameter<std::string>("object");
0081   labelName_ = iConfig.getParameter<std::string>("label");
0082   if (objectType_ == "electron")
0083     electronsToken_ = consumes<std::vector<pat::Electron> >(edm::InputTag(labelName_));
0084   else if (objectType_ == "muon")
0085     muonsToken_ = consumes<std::vector<pat::Muon> >(edm::InputTag(labelName_));
0086   else if (objectType_ == "lJets" || objectType_ == "bJets")
0087     jetsToken_ = consumes<std::vector<pat::Jet> >(edm::InputTag(labelName_));
0088   else if (objectType_ == "met")
0089     metsToken_ = consumes<std::vector<pat::MET> >(edm::InputTag(labelName_));
0090   else if (objectType_ == "tau")
0091     tausToken_ = consumes<std::vector<pat::Tau> >(edm::InputTag(labelName_));
0092   if (objectType_ != "met") {
0093     etabinVals_ = iConfig.getParameter<std::vector<double> >("etabinValues");
0094   }
0095   pTbinVals_ = iConfig.getParameter<std::vector<double> >("pTbinValues");
0096   minDR_ = iConfig.getParameter<double>("minMatchingDR");
0097 
0098   nrFilled = 0;
0099 }
0100 
0101 //
0102 // member functions
0103 //
0104 
0105 // ------------ method called to for each event  ------------
0106 void ResolutionCreator::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0107   // Get the gen and cal object fourvector
0108   std::vector<reco::Particle *> p4gen, p4rec;
0109 
0110   edm::Handle<TtGenEvent> genEvt;
0111   iEvent.getByToken(genEvtToken_, genEvt);
0112 
0113   if (genEvt->particles().size() < 10)
0114     return;
0115 
0116   if (objectType_ == "electron") {
0117     edm::Handle<std::vector<pat::Electron> >
0118         electrons;  //to calculate the ResolutionCreator for the electrons, i used the TopElectron instead of the AOD information
0119     iEvent.getByToken(electronsToken_, electrons);
0120     for (size_t e = 0; e < electrons->size(); e++) {
0121       for (size_t p = 0; p < genEvt->particles().size(); p++) {
0122         if ((std::abs(genEvt->particles()[p].pdgId()) == 11) &&
0123             (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*electrons)[e].p4()) < minDR_)) {
0124           //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
0125           //p4rec.push_back(new reco::Particle((pat::Electron)((*electrons)[e])));
0126         }
0127       }
0128     }
0129   } else if (objectType_ == "muon") {
0130     edm::Handle<std::vector<pat::Muon> > muons;
0131     iEvent.getByToken(muonsToken_, muons);
0132     for (size_t m = 0; m < muons->size(); m++) {
0133       for (size_t p = 0; p < genEvt->particles().size(); p++) {
0134         if ((std::abs(genEvt->particles()[p].pdgId()) == 13) &&
0135             (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*muons)[m].p4()) < minDR_)) {
0136           //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
0137           //p4rec.push_back(new reco::Particle((pat::Muon)((*muons)[m])));
0138         }
0139       }
0140     }
0141   } else if (objectType_ == "lJets") {
0142     edm::Handle<std::vector<pat::Jet> > jets;
0143     iEvent.getByToken(jetsToken_, jets);
0144     if (jets->size() >= 4) {
0145       for (unsigned int j = 0; j < 4; j++) {
0146         for (size_t p = 0; p < genEvt->particles().size(); p++) {
0147           if ((std::abs(genEvt->particles()[p].pdgId()) < 5) &&
0148               (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*jets)[j].p4()) < minDR_)) {
0149             //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
0150             //p4rec.push_back(new reco::Particle((pat::Jet)(*jets)[j]));
0151           }
0152         }
0153       }
0154     }
0155   } else if (objectType_ == "bJets") {
0156     edm::Handle<std::vector<pat::Jet> > jets;
0157     iEvent.getByToken(jetsToken_, jets);
0158     if (jets->size() >= 4) {
0159       for (unsigned int j = 0; j < 4; j++) {
0160         for (size_t p = 0; p < genEvt->particles().size(); p++) {
0161           if ((std::abs(genEvt->particles()[p].pdgId()) == 5) &&
0162               (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*jets)[j].p4()) < minDR_)) {
0163             //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
0164             //p4rec.push_back(new reco::Particle((pat::Jet)(*jets)[j]));
0165           }
0166         }
0167       }
0168     }
0169   } else if (objectType_ == "met") {
0170     edm::Handle<std::vector<pat::MET> > mets;
0171     iEvent.getByToken(metsToken_, mets);
0172     if (!mets->empty()) {
0173       if (genEvt->isSemiLeptonic() && genEvt->singleNeutrino() != nullptr &&
0174           ROOT::Math::VectorUtil::DeltaR(genEvt->singleNeutrino()->p4(), (*mets)[0].p4()) < minDR_) {
0175         //p4gen.push_back(new reco::Particle(0,genEvt->singleNeutrino()->p4(),math::XYZPoint()));
0176         //p4rec.push_back(new reco::Particle((pat::MET)((*mets)[0])));
0177       }
0178     }
0179   } else if (objectType_ == "tau") {
0180     edm::Handle<std::vector<pat::Tau> > taus;
0181     iEvent.getByToken(tausToken_, taus);
0182     for (std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau != taus->end(); ++tau) {
0183       // find the tau (if any) that matches a MC tau from W
0184       reco::GenParticle genLepton = *(tau->genLepton());
0185       if (std::abs(genLepton.pdgId()) == 15 && genLepton.status() == 2 && genLepton.numberOfMothers() > 0 &&
0186           std::abs(genLepton.mother(0)->pdgId()) == 15 && genLepton.mother(0)->numberOfMothers() > 0 &&
0187           std::abs(genLepton.mother(0)->mother(0)->pdgId()) == 24 &&
0188           ROOT::Math::VectorUtil::DeltaR(genLepton.p4(), tau->p4()) < minDR_) {
0189       }
0190       //p4gen.push_back(new reco::Particle(genLepton));
0191       //p4rec.push_back(new reco::Particle(*tau));
0192     }
0193   }
0194   // Fill the object's value
0195   for (unsigned m = 0; m < p4gen.size(); m++) {
0196     double Egen = p4gen[m]->energy();
0197     double Thetagen = p4gen[m]->theta();
0198     double Phigen = p4gen[m]->phi();
0199     double Etgen = p4gen[m]->et();
0200     double Etagen = p4gen[m]->eta();
0201     double Ecal = p4rec[m]->energy();
0202     double Thetacal = p4rec[m]->theta();
0203     double Phical = p4rec[m]->phi();
0204     double Etcal = p4rec[m]->et();
0205     double Etacal = p4rec[m]->eta();
0206     double phidiff = Phical - Phigen;
0207     if (phidiff > 3.14159)
0208       phidiff = 2. * 3.14159 - phidiff;
0209     if (phidiff < -3.14159)
0210       phidiff = -phidiff - 2. * 3.14159;
0211 
0212     // find eta and et bin
0213     int etabin = 0;
0214     if (etanrbins > 1) {
0215       for (unsigned int b = 0; b < etabinVals_.size() - 1; b++) {
0216         if (fabs(Etacal) > etabinVals_[b])
0217           etabin = b;
0218       }
0219     }
0220 
0221     int ptbin = 0;
0222     for (unsigned int b = 0; b < pTbinVals_.size() - 1; b++) {
0223       if (p4rec[m]->pt() > pTbinVals_[b])
0224         ptbin = b;
0225     }
0226 
0227     // calculate the resolution on "a", "b", "c" & "d" according to the definition (CMS-NOTE-2006-023):
0228     // p = a*|p_meas|*u_1 + b*u_2 + c*u_3
0229     // E(fit) = E_meas * d
0230     //
0231     // with u_1 = p/|p_meas|
0232     //      u_3 = (u_z x u_1)/|u_z x u_1|
0233     //      u_2 = (u_1 x u_3)/|u_1 x u_3|
0234     //
0235     // The initial parameters values are chosen like (a, b, c, d) = (1., 0., 0., 1.)
0236 
0237     // 1/ calculate the unitary vectors of the basis u_1, u_2, u_3
0238     ROOT::Math::SVector<double, 3> pcalvec(p4rec[m]->px(), p4rec[m]->py(), p4rec[m]->pz());
0239     ROOT::Math::SVector<double, 3> pgenvec(p4gen[m]->px(), p4gen[m]->py(), p4gen[m]->pz());
0240 
0241     ROOT::Math::SVector<double, 3> u_z(0, 0, 1);
0242     ROOT::Math::SVector<double, 3> u_1 = ROOT::Math::Unit(pcalvec);
0243     ROOT::Math::SVector<double, 3> u_3 = ROOT::Math::Cross(u_z, u_1) / ROOT::Math::Mag(ROOT::Math::Cross(u_z, u_1));
0244     ROOT::Math::SVector<double, 3> u_2 = ROOT::Math::Cross(u_1, u_3) / ROOT::Math::Mag(ROOT::Math::Cross(u_1, u_3));
0245     double acal = 1.;
0246     double bcal = 0.;
0247     double ccal = 0.;
0248     double dcal = 1.;
0249     double agen = ROOT::Math::Dot(pgenvec, u_1) / ROOT::Math::Mag(pcalvec);
0250     double bgen = ROOT::Math::Dot(pgenvec, u_2);
0251     double cgen = ROOT::Math::Dot(pgenvec, u_3);
0252     double dgen = Egen / Ecal;
0253 
0254     //fill histograms
0255     ++nrFilled;
0256     hResPtEtaBin[0][etabin][ptbin]->Fill(acal - agen);
0257     hResPtEtaBin[1][etabin][ptbin]->Fill(bcal - bgen);
0258     hResPtEtaBin[2][etabin][ptbin]->Fill(ccal - cgen);
0259     hResPtEtaBin[3][etabin][ptbin]->Fill(dcal - dgen);
0260     hResPtEtaBin[4][etabin][ptbin]->Fill(Thetacal - Thetagen);
0261     hResPtEtaBin[5][etabin][ptbin]->Fill(phidiff);
0262     hResPtEtaBin[6][etabin][ptbin]->Fill(Etcal - Etgen);
0263     hResPtEtaBin[7][etabin][ptbin]->Fill(Etacal - Etagen);
0264 
0265     delete p4gen[m];
0266     delete p4rec[m];
0267   }
0268 }
0269 
0270 // ------------ method called once each job just before starting event loop  ------------
0271 void ResolutionCreator::beginJob() {
0272   edm::Service<TFileService> fs;
0273   if (!fs)
0274     throw edm::Exception(edm::errors::Configuration, "TFileService missing from configuration!");
0275 
0276   // input constants
0277   TString resObsName[8] = {"_ares", "_bres", "_cres", "_dres", "_thres", "_phres", "_etres", "_etares"};
0278   int resObsNrBins = 120;
0279   if ((objectType_ == "muon") || (objectType_ == "electron"))
0280     resObsNrBins = 80;
0281   std::vector<double> resObsMin, resObsMax;
0282   if (objectType_ == "electron") {
0283     resObsMin.push_back(-0.15);
0284     resObsMin.push_back(-0.2);
0285     resObsMin.push_back(-0.1);
0286     resObsMin.push_back(-0.15);
0287     resObsMin.push_back(-0.0012);
0288     resObsMin.push_back(-0.009);
0289     resObsMin.push_back(-16);
0290     resObsMin.push_back(-0.0012);
0291     resObsMax.push_back(0.15);
0292     resObsMax.push_back(0.2);
0293     resObsMax.push_back(0.1);
0294     resObsMax.push_back(0.15);
0295     resObsMax.push_back(0.0012);
0296     resObsMax.push_back(0.009);
0297     resObsMax.push_back(16);
0298     resObsMax.push_back(0.0012);
0299   } else if (objectType_ == "muon") {
0300     resObsMin.push_back(-0.15);
0301     resObsMin.push_back(-0.1);
0302     resObsMin.push_back(-0.05);
0303     resObsMin.push_back(-0.15);
0304     resObsMin.push_back(-0.004);
0305     resObsMin.push_back(-0.003);
0306     resObsMin.push_back(-8);
0307     resObsMin.push_back(-0.004);
0308     resObsMax.push_back(0.15);
0309     resObsMax.push_back(0.1);
0310     resObsMax.push_back(0.05);
0311     resObsMax.push_back(0.15);
0312     resObsMax.push_back(0.004);
0313     resObsMax.push_back(0.003);
0314     resObsMax.push_back(8);
0315     resObsMax.push_back(0.004);
0316   } else if (objectType_ == "tau") {
0317     resObsMin.push_back(-1.);
0318     resObsMin.push_back(-10.);
0319     resObsMin.push_back(-10);
0320     resObsMin.push_back(-1.);
0321     resObsMin.push_back(-0.1);
0322     resObsMin.push_back(-0.1);
0323     resObsMin.push_back(-80);
0324     resObsMin.push_back(-0.1);
0325     resObsMax.push_back(1.);
0326     resObsMax.push_back(10.);
0327     resObsMax.push_back(10);
0328     resObsMax.push_back(1.);
0329     resObsMax.push_back(0.1);
0330     resObsMax.push_back(0.1);
0331     resObsMax.push_back(50);
0332     resObsMax.push_back(0.1);
0333   } else if (objectType_ == "lJets" || objectType_ == "bJets") {
0334     resObsMin.push_back(-1.);
0335     resObsMin.push_back(-10.);
0336     resObsMin.push_back(-10.);
0337     resObsMin.push_back(-1.);
0338     resObsMin.push_back(-0.4);
0339     resObsMin.push_back(-0.6);
0340     resObsMin.push_back(-80);
0341     resObsMin.push_back(-0.6);
0342     resObsMax.push_back(1.);
0343     resObsMax.push_back(10.);
0344     resObsMax.push_back(10.);
0345     resObsMax.push_back(1.);
0346     resObsMax.push_back(0.4);
0347     resObsMax.push_back(0.6);
0348     resObsMax.push_back(80);
0349     resObsMax.push_back(0.6);
0350   } else {
0351     resObsMin.push_back(-2.);
0352     resObsMin.push_back(-150.);
0353     resObsMin.push_back(-150.);
0354     resObsMin.push_back(-2.);
0355     resObsMin.push_back(-6);
0356     resObsMin.push_back(-6);
0357     resObsMin.push_back(-180);
0358     resObsMin.push_back(-6);
0359     resObsMax.push_back(3.);
0360     resObsMax.push_back(150.);
0361     resObsMax.push_back(150.);
0362     resObsMax.push_back(3.);
0363     resObsMax.push_back(6);
0364     resObsMax.push_back(6);
0365     resObsMax.push_back(180);
0366     resObsMax.push_back(6);
0367   }
0368 
0369   const char *resObsVsPtFit[8] = {"[0]+[1]*exp(-[2]*x)",
0370                                   "[0]+[1]*exp(-[2]*x)",
0371                                   "[0]+[1]*exp(-[2]*x)",
0372                                   "[0]+[1]*exp(-[2]*x)",
0373                                   "[0]+[1]*exp(-[2]*x)",
0374                                   "[0]+[1]*exp(-[2]*x)",
0375                                   "pol1",
0376                                   "[0]+[1]*exp(-[2]*x)"};
0377 
0378   ptnrbins = pTbinVals_.size() - 1;
0379   double *ptbins = new double[pTbinVals_.size()];
0380   for (unsigned int b = 0; b < pTbinVals_.size(); b++)
0381     ptbins[b] = pTbinVals_[b];
0382   double *etabins;
0383   if (objectType_ != "met") {
0384     etanrbins = etabinVals_.size() - 1;
0385     etabins = new double[etabinVals_.size()];
0386     for (unsigned int b = 0; b < etabinVals_.size(); b++)
0387       etabins[b] = etabinVals_[b];
0388   } else {
0389     etanrbins = 1;
0390     etabins = new double[2];
0391     etabins[0] = 0;
0392     etabins[1] = 5.;
0393   }
0394 
0395   //define the histograms booked
0396   for (Int_t ro = 0; ro < 8; ro++) {
0397     for (Int_t etab = 0; etab < etanrbins; etab++) {
0398       for (Int_t ptb = 0; ptb < ptnrbins; ptb++) {
0399         TString obsName = objectType_;
0400         obsName += resObsName[ro];
0401         obsName += "_etabin";
0402         obsName += etab;
0403         obsName += "_ptbin";
0404         obsName += ptb;
0405         hResPtEtaBin[ro][etab][ptb] = fs->make<TH1F>(obsName, obsName, resObsNrBins, resObsMin[ro], resObsMax[ro]);
0406         fResPtEtaBin[ro][etab][ptb] = fs->make<TF1>("F_" + obsName, "gaus");
0407       }
0408       TString obsName2 = objectType_;
0409       obsName2 += resObsName[ro];
0410       obsName2 += "_etabin";
0411       obsName2 += etab;
0412       hResEtaBin[ro][etab] = fs->make<TH1F>(obsName2, obsName2, ptnrbins, ptbins);
0413       fResEtaBin[ro][etab] =
0414           fs->make<TF1>("F_" + obsName2, resObsVsPtFit[ro], pTbinVals_[0], pTbinVals_[pTbinVals_.size() - 1]);
0415     }
0416   }
0417   tResVar = fs->make<TTree>("tResVar", "Resolution tree");
0418 
0419   delete[] etabins;
0420   delete[] ptbins;
0421 }
0422 
0423 // ------------ method called once each job just after ending the event loop  ------------
0424 void ResolutionCreator::endJob() {
0425   TString resObsName2[8] = {"a", "b", "c", "d", "theta", "phi", "et", "eta"};
0426   Int_t ro = 0;
0427   Double_t pt = 0.;
0428   Double_t eta = 0.;
0429   Double_t value, error;
0430 
0431   tResVar->Branch("Pt", &pt, "Pt/D");
0432   tResVar->Branch("Eta", &eta, "Eta/D");
0433   tResVar->Branch("ro", &ro, "ro/I");
0434   tResVar->Branch("value", &value, "value/D");
0435   tResVar->Branch("error", &error, "error/D");
0436 
0437   for (ro = 0; ro < 8; ro++) {
0438     for (int etab = 0; etab < etanrbins; etab++) {
0439       //CD set eta at the center of the bin
0440       eta = etanrbins > 1 ? (etabinVals_[etab] + etabinVals_[etab + 1]) / 2. : 2.5;
0441       for (int ptb = 0; ptb < ptnrbins; ptb++) {
0442         //CD set pt at the center of the bin
0443         pt = (pTbinVals_[ptb] + pTbinVals_[ptb + 1]) / 2.;
0444         double maxcontent = 0.;
0445         int maxbin = 0;
0446         for (int nb = 1; nb < hResPtEtaBin[ro][etab][ptb]->GetNbinsX(); nb++) {
0447           if (hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb) > maxcontent) {
0448             maxcontent = hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb);
0449             maxbin = nb;
0450           }
0451         }
0452         int range = (int)(hResPtEtaBin[ro][etab][ptb]->GetNbinsX() / 6);  //in order that ~1/3 of X-axis range is fitted
0453         fResPtEtaBin[ro][etab][ptb]->SetRange(hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin - range),
0454                                               hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin + range));
0455         fResPtEtaBin[ro][etab][ptb]->SetParameters(hResPtEtaBin[ro][etab][ptb]->GetMaximum(),
0456                                                    hResPtEtaBin[ro][etab][ptb]->GetMean(),
0457                                                    hResPtEtaBin[ro][etab][ptb]->GetRMS());
0458         hResPtEtaBin[ro][etab][ptb]->Fit(fResPtEtaBin[ro][etab][ptb]->GetName(), "RQ SERIAL");
0459         hResEtaBin[ro][etab]->SetBinContent(ptb + 1, fResPtEtaBin[ro][etab][ptb]->GetParameter(2));
0460         hResEtaBin[ro][etab]->SetBinError(ptb + 1, fResPtEtaBin[ro][etab][ptb]->GetParError(2));
0461         //CD: Fill the tree
0462         value = fResPtEtaBin[ro][etab][ptb]->GetParameter(2);  //parameter value
0463         error = fResPtEtaBin[ro][etab][ptb]->GetParError(2);   //parameter error
0464         tResVar->Fill();
0465       }
0466       //CD: add a fake entry in pt=0 for the NN training
0467       // for that, use a linear extrapolation.
0468       pt = 0.;
0469       value = ((pTbinVals_[0] + pTbinVals_[1]) / 2.) *
0470                   (fResPtEtaBin[ro][etab][0]->GetParameter(2) - fResPtEtaBin[ro][etab][1]->GetParameter(2)) /
0471                   ((pTbinVals_[2] - pTbinVals_[0]) / 2.) +
0472               fResPtEtaBin[ro][etab][0]->GetParameter(2);
0473       error = fResPtEtaBin[ro][etab][0]->GetParError(2) + fResPtEtaBin[ro][etab][1]->GetParError(2);
0474       tResVar->Fill();
0475       // standard fit
0476       hResEtaBin[ro][etab]->Fit(fResEtaBin[ro][etab]->GetName(), "RQ SERIAL");
0477     }
0478   }
0479   if (objectType_ == "lJets" && nrFilled == 0)
0480     edm::LogProblem("SummaryError") << "No plots filled for light jets \n";
0481   if (objectType_ == "bJets" && nrFilled == 0)
0482     edm::LogProblem("SummaryError") << "No plots filled for bjets \n";
0483   if (objectType_ == "muon" && nrFilled == 0)
0484     edm::LogProblem("SummaryError") << "No plots filled for muons \n";
0485   if (objectType_ == "electron" && nrFilled == 0)
0486     edm::LogProblem("SummaryError") << "No plots filled for electrons \n";
0487   if (objectType_ == "tau" && nrFilled == 0)
0488     edm::LogProblem("SummaryError") << "No plots filled for taus \n";
0489   if (objectType_ == "met" && nrFilled == 0)
0490     edm::LogProblem("SummaryError") << "No plots filled for met \n";
0491 
0492   edm::LogVerbatim("MainResults") << " \n\n";
0493   edm::LogVerbatim("MainResults") << " ----------------------------------------------";
0494   edm::LogVerbatim("MainResults") << " ----------------------------------------------";
0495   edm::LogVerbatim("MainResults") << " Resolutions on " << objectType_ << " with nrfilled: " << nrFilled;
0496   edm::LogVerbatim("MainResults") << " ----------------------------------------------";
0497   edm::LogVerbatim("MainResults") << " ----------------------------------------------";
0498   if (nrFilled != 0 && objectType_ != "met") {
0499     for (ro = 0; ro < 8; ro++) {
0500       edm::LogVerbatim("MainResults") << "-------------------- ";
0501       edm::LogVerbatim("MainResults") << "\n Resolutions on " << resObsName2[ro] << "\n";
0502       edm::LogVerbatim("MainResults") << "-------------------- ";
0503       for (int etab = 0; etab < etanrbins; etab++) {
0504         if (nrFilled != 0 && ro != 6) {
0505           if (etab == 0) {
0506             edm::LogVerbatim("MainResults")
0507                 << "if(fabs(eta)<" << etabinVals_[etab + 1] << ") res = " << fResEtaBin[ro][etab]->GetParameter(0)
0508                 << "+" << fResEtaBin[ro][etab]->GetParameter(1) << "*exp(-(" << fResEtaBin[ro][etab]->GetParameter(2)
0509                 << "*pt));";
0510           } else {
0511             edm::LogVerbatim("MainResults")
0512                 << "else if(fabs(eta)<" << etabinVals_[etab + 1] << ") res = " << fResEtaBin[ro][etab]->GetParameter(0)
0513                 << "+" << fResEtaBin[ro][etab]->GetParameter(1) << "*exp(-(" << fResEtaBin[ro][etab]->GetParameter(2)
0514                 << "*pt));";
0515           }
0516         } else if (nrFilled != 0 && ro == 6) {
0517           if (etab == 0) {
0518             edm::LogVerbatim("MainResults")
0519                 << "if(fabs(eta)<" << etabinVals_[etab + 1] << ") res = " << fResEtaBin[ro][etab]->GetParameter(0)
0520                 << "+" << fResEtaBin[ro][etab]->GetParameter(1) << "*pt;";
0521           } else {
0522             edm::LogVerbatim("MainResults")
0523                 << "else if(fabs(eta)<" << etabinVals_[etab + 1] << ") res = " << fResEtaBin[ro][etab]->GetParameter(0)
0524                 << "+" << fResEtaBin[ro][etab]->GetParameter(1) << "*pt;";
0525           }
0526         }
0527       }
0528     }
0529   } else if (nrFilled != 0 && objectType_ == "met") {
0530     for (ro = 0; ro < 8; ro++) {
0531       edm::LogVerbatim("MainResults") << "-------------------- ";
0532       edm::LogVerbatim("MainResults") << "\n Resolutions on " << resObsName2[ro] << "\n";
0533       edm::LogVerbatim("MainResults") << "-------------------- ";
0534       if (nrFilled != 0 && ro != 6) {
0535         edm::LogVerbatim("MainResults") << "res = " << fResEtaBin[ro][0]->GetParameter(0) << "+"
0536                                         << fResEtaBin[ro][0]->GetParameter(1) << "*exp(-("
0537                                         << fResEtaBin[ro][0]->GetParameter(2) << "*pt));";
0538       } else if (nrFilled != 0 && ro == 6) {
0539         edm::LogVerbatim("MainResults") << "res = " << fResEtaBin[ro][0]->GetParameter(0) << "+"
0540                                         << fResEtaBin[ro][0]->GetParameter(1) << "*pt;";
0541       }
0542     }
0543   }
0544 }
0545 
0546 //define this as a plug-in
0547 DEFINE_FWK_MODULE(ResolutionCreator);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!