CMS 3D CMS Logo
ResolutionCreator.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 #include <string>
4 
5 // user include files
6 #include "FWCore/Framework/interface/Frameworkfwd.h"
7 #include "FWCore/Framework/interface/EDAnalyzer.h"
8 
9 #include "FWCore/Framework/interface/Event.h"
10 #include "FWCore/Framework/interface/MakerMacros.h"
11 
12 #include "FWCore/ParameterSet/interface/ParameterSet.h"
13 //needed for TFileService
14 #include "CommonTools/UtilAlgos/interface/TFileService.h"
15 #include "FWCore/ServiceRegistry/interface/Service.h"
16 //needed for MessageLogger
17 #include "FWCore/MessageLogger/interface/MessageLogger.h"
18 
19 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
20 
21 #include "DataFormats/PatCandidates/interface/Muon.h"
22 #include "DataFormats/PatCandidates/interface/Jet.h"
23 #include "DataFormats/PatCandidates/interface/MET.h"
24 #include "DataFormats/PatCandidates/interface/Tau.h"
25 #include "DataFormats/PatCandidates/interface/Electron.h"
26 
27 #include "TDirectory.h"
28 #include "TH1F.h"
29 #include "TF1.h"
30 #include "TTree.h"
31 
32 #include <Math/VectorUtil.h>
33 
34 //
35 // class declaration
36 //
37 
38 class ResolutionCreator : public edm::EDAnalyzer {
39  public:
40  explicit ResolutionCreator(const edm::ParameterSet&);
41  ~ResolutionCreator() override;
42 
43  private:
44  void beginJob() override ;
45  void analyze(const edm::Event&, const edm::EventSetup&) override;
46  void endJob() override ;
47 
48  // ----------member data ---------------------------
49  edm::EDGetTokenT<TtGenEvent> genEvtToken_;
50  std::string objectType_, labelName_;
51  edm::EDGetTokenT<std::vector<pat::Electron> > electronsToken_;
52  edm::EDGetTokenT<std::vector<pat::Muon> > muonsToken_;
53  edm::EDGetTokenT<std::vector<pat::Jet> > jetsToken_;
54  edm::EDGetTokenT<std::vector<pat::MET> > metsToken_;
55  edm::EDGetTokenT<std::vector<pat::Tau> > tausToken_;
56  std::vector<double> etabinVals_, pTbinVals_;
57  double minDR_;
58  int matchingAlgo_;
59  bool useMaxDist_;
60  bool useDeltaR_;
61  double maxDist_;
62  int ptnrbins, etanrbins;
63  int nrFilled;
64 
65  //Histograms are booked in the beginJob() method
66  TF1 *fResPtEtaBin[10][20][20];
67  TF1 *fResEtaBin[10][20];
68  TH1F *hResPtEtaBin[10][20][20];
69  TH1F *hResEtaBin[10][20];
70  TTree* tResVar;
71 };
72 
73 
74 //
75 // constructors and destructor
76 //
77 ResolutionCreator::ResolutionCreator(const edm::ParameterSet& iConfig)
78 {
79  // input parameters
80  genEvtToken_ = consumes<TtGenEvent>(edm::InputTag("genEvt"));
81  objectType_ = iConfig.getParameter< std::string > ("object");
82  labelName_ = iConfig.getParameter< std::string > ("label");
83  if(objectType_ == "electron") electronsToken_ = consumes<std::vector<pat::Electron> >(edm::InputTag(labelName_));
84  else if(objectType_ == "muon") muonsToken_ = consumes<std::vector<pat::Muon> >(edm::InputTag(labelName_));
85  else if(objectType_ == "lJets" || objectType_ == "bJets") jetsToken_ = consumes<std::vector<pat::Jet> >(edm::InputTag(labelName_));
86  else if(objectType_ == "met") metsToken_ = consumes<std::vector<pat::MET> >(edm::InputTag(labelName_));
87  else if(objectType_ == "tau") tausToken_ = consumes<std::vector<pat::Tau> >(edm::InputTag(labelName_));
88  if(objectType_ != "met"){
89  etabinVals_ = iConfig.getParameter< std::vector<double> > ("etabinValues");
90  }
91  pTbinVals_ = iConfig.getParameter< std::vector<double> > ("pTbinValues");
92  minDR_ = iConfig.getParameter< double > ("minMatchingDR");
93 
94  nrFilled = 0;
95 
96 }
97 
98 
99 ResolutionCreator::~ResolutionCreator()
100 {
101 }
102 
103 
104 //
105 // member functions
106 //
107 
108 // ------------ method called to for each event ------------
109 void
110 ResolutionCreator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
111 {
112 
113  // Get the gen and cal object fourvector
114  std::vector<reco::Particle *> p4gen, p4rec;
115 
116  edm::Handle<TtGenEvent> genEvt;
117  iEvent.getByToken(genEvtToken_,genEvt);
118 
119  if(genEvt->particles().size()<10) return;
120 
121  if(objectType_ == "electron"){
122  edm::Handle<std::vector<pat::Electron> > electrons; //to calculate the ResolutionCreator for the electrons, i used the TopElectron instead of the AOD information
123  iEvent.getByToken(electronsToken_,electrons);
124  for(size_t e=0; e<electrons->size(); e++) {
125  for(size_t p=0; p<genEvt->particles().size(); p++){
126  if( (std::abs(genEvt->particles()[p].pdgId()) == 11) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*electrons)[e].p4()) < minDR_) ) {
127  //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
128  //p4rec.push_back(new reco::Particle((pat::Electron)((*electrons)[e])));
129  }
130  }
131  }
132  }
133  else if(objectType_ == "muon"){
134  edm::Handle<std::vector<pat::Muon> > muons;
135  iEvent.getByToken(muonsToken_,muons);
136  for(size_t m=0; m<muons->size(); m++) {
137  for(size_t p=0; p<genEvt->particles().size(); p++){
138  if( (std::abs(genEvt->particles()[p].pdgId()) == 13) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*muons)[m].p4()) < minDR_) ) {
139  //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
140  //p4rec.push_back(new reco::Particle((pat::Muon)((*muons)[m])));
141  }
142  }
143  }
144  }
145  else if(objectType_ == "lJets" ){
146  edm::Handle<std::vector<pat::Jet> > jets;
147  iEvent.getByToken(jetsToken_,jets);
148  if(jets->size()>=4) {
149  for(unsigned int j = 0; j<4; j++){
150  for(size_t p=0; p<genEvt->particles().size(); p++){
151  if( (std::abs(genEvt->particles()[p].pdgId()) < 5) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*jets)[j].p4())< minDR_) ){
152  //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
153  //p4rec.push_back(new reco::Particle((pat::Jet)(*jets)[j]));
154  }
155  }
156  }
157  }
158  }
159  else if(objectType_ == "bJets" ){
160  edm::Handle<std::vector<pat::Jet> > jets;
161  iEvent.getByToken(jetsToken_,jets);
162  if(jets->size()>=4) {
163  for(unsigned int j = 0; j<4; j++){
164  for(size_t p=0; p<genEvt->particles().size(); p++){
165  if( (std::abs(genEvt->particles()[p].pdgId()) == 5) && (ROOT::Math::VectorUtil::DeltaR(genEvt->particles()[p].p4(), (*jets)[j].p4())< minDR_) ) {
166  //p4gen.push_back(new reco::Particle(genEvt->particles()[p]));
167  //p4rec.push_back(new reco::Particle((pat::Jet)(*jets)[j]));
168  }
169  }
170  }
171  }
172  }
173  else if(objectType_ == "met"){
174  edm::Handle<std::vector<pat::MET> > mets;
175  iEvent.getByToken(metsToken_,mets);
176  if(!mets->empty()) {
177  if( genEvt->isSemiLeptonic() && genEvt->singleNeutrino() != nullptr && ROOT::Math::VectorUtil::DeltaR(genEvt->singleNeutrino()->p4(), (*mets)[0].p4()) < minDR_) {
178  //p4gen.push_back(new reco::Particle(0,genEvt->singleNeutrino()->p4(),math::XYZPoint()));
179  //p4rec.push_back(new reco::Particle((pat::MET)((*mets)[0])));
180  }
181  }
182  }
183  else if(objectType_ == "tau"){
184  edm::Handle<std::vector<pat::Tau> > taus;
185  iEvent.getByToken(tausToken_,taus);
186  for(std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau != taus->end(); ++tau) {
187  // find the tau (if any) that matches a MC tau from W
188  reco::GenParticle genLepton = *(tau->genLepton());
189  if( std::abs(genLepton.pdgId())==15 && genLepton.status()==2 &&
190  genLepton.numberOfMothers()>0 &&
191  std::abs(genLepton.mother(0)->pdgId())==15 &&
192  genLepton.mother(0)->numberOfMothers()>0 &&
193  std::abs(genLepton.mother(0)->mother(0)->pdgId())==24 &&
194  ROOT::Math::VectorUtil::DeltaR(genLepton.p4(), tau->p4()) < minDR_ ) {
195  }
196  //p4gen.push_back(new reco::Particle(genLepton));
197  //p4rec.push_back(new reco::Particle(*tau));
198  }
199  }
200  // Fill the object's value
201  for(unsigned m=0; m<p4gen.size(); m++){
202  double Egen = p4gen[m]->energy();
203  double Thetagen = p4gen[m]->theta();
204  double Phigen = p4gen[m]->phi();
205  double Etgen = p4gen[m]->et();
206  double Etagen = p4gen[m]->eta();
207  double Ecal = p4rec[m]->energy();
208  double Thetacal = p4rec[m]->theta();
209  double Phical = p4rec[m]->phi();
210  double Etcal = p4rec[m]->et();
211  double Etacal = p4rec[m]->eta();
212  double phidiff = Phical- Phigen;
213  if(phidiff>3.14159) phidiff = 2.*3.14159 - phidiff;
214  if(phidiff<-3.14159) phidiff = -phidiff - 2.*3.14159;
215 
216  // find eta and et bin
217  int etabin = 0;
218  if(etanrbins > 1){
219  for(unsigned int b=0; b<etabinVals_.size()-1; b++) {
220  if(fabs(Etacal) > etabinVals_[b]) etabin = b;
221  }
222  }
223 
224  int ptbin = 0;
225  for(unsigned int b=0; b<pTbinVals_.size()-1; b++) {
226  if(p4rec[m]->pt() > pTbinVals_[b]) ptbin = b;
227  }
228 
229  // calculate the resolution on "a", "b", "c" & "d" according to the definition (CMS-NOTE-2006-023):
230  // p = a*|p_meas|*u_1 + b*u_2 + c*u_3
231  // E(fit) = E_meas * d
232  //
233  // with u_1 = p/|p_meas|
234  // u_3 = (u_z x u_1)/|u_z x u_1|
235  // u_2 = (u_1 x u_3)/|u_1 x u_3|
236  //
237  // The initial parameters values are chosen like (a, b, c, d) = (1., 0., 0., 1.)
238 
239  // 1/ calculate the unitary vectors of the basis u_1, u_2, u_3
240  ROOT::Math::SVector<double,3> pcalvec(p4rec[m]->px(),p4rec[m]->py(),p4rec[m]->pz());
241  ROOT::Math::SVector<double,3> pgenvec(p4gen[m]->px(),p4gen[m]->py(),p4gen[m]->pz());
242 
243  ROOT::Math::SVector<double,3> u_z(0,0,1);
244  ROOT::Math::SVector<double,3> u_1 = ROOT::Math::Unit(pcalvec);
245  ROOT::Math::SVector<double,3> u_3 = ROOT::Math::Cross(u_z,u_1)/ROOT::Math::Mag(ROOT::Math::Cross(u_z,u_1));
246  ROOT::Math::SVector<double,3> u_2 = ROOT::Math::Cross(u_1,u_3)/ROOT::Math::Mag(ROOT::Math::Cross(u_1,u_3));
247  double acal = 1.;
248  double bcal = 0.;
249  double ccal = 0.;
250  double dcal = 1.;
251  double agen = ROOT::Math::Dot(pgenvec,u_1)/ROOT::Math::Mag(pcalvec);
252  double bgen = ROOT::Math::Dot(pgenvec,u_2);
253  double cgen = ROOT::Math::Dot(pgenvec,u_3);
254  double dgen = Egen/Ecal;
255 
256  //fill histograms
257  ++nrFilled;
258  hResPtEtaBin[0][etabin][ptbin] -> Fill(acal-agen);
259  hResPtEtaBin[1][etabin][ptbin] -> Fill(bcal-bgen);
260  hResPtEtaBin[2][etabin][ptbin] -> Fill(ccal-cgen);
261  hResPtEtaBin[3][etabin][ptbin] -> Fill(dcal-dgen);
262  hResPtEtaBin[4][etabin][ptbin] -> Fill(Thetacal-Thetagen);
263  hResPtEtaBin[5][etabin][ptbin] -> Fill(phidiff);
264  hResPtEtaBin[6][etabin][ptbin] -> Fill(Etcal-Etgen);
265  hResPtEtaBin[7][etabin][ptbin] -> Fill(Etacal-Etagen);
266 
267  delete p4gen[m];
268  delete p4rec[m];
269  }
270 
271 }
272 
273 
274 // ------------ method called once each job just before starting event loop ------------
275 void
276 ResolutionCreator::beginJob()
277 {
278  edm::Service<TFileService> fs;
279  if (!fs) throw edm::Exception(edm::errors::Configuration, "TFileService missing from configuration!");
280 
281  // input constants
282  TString resObsName[8] = {"_ares","_bres","_cres","_dres","_thres","_phres","_etres","_etares"};
283  int resObsNrBins = 120;
284  if( (objectType_ == "muon") || (objectType_ == "electron") ) resObsNrBins = 80;
285  std::vector<double> resObsMin, resObsMax;
286  if(objectType_ == "electron"){
287  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);
288  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);
289  } else if(objectType_ == "muon"){
290  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);
291  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);
292  } else if(objectType_ == "tau"){
293  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);
294  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);
295  } else if(objectType_ == "lJets" || objectType_ == "bJets"){
296  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);
297  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);
298  } else{
299  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);
300  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);
301  }
302 
303  const char* resObsVsPtFit[8] = { "[0]+[1]*exp(-[2]*x)",
304  "[0]+[1]*exp(-[2]*x)",
305  "[0]+[1]*exp(-[2]*x)",
306  "[0]+[1]*exp(-[2]*x)",
307  "[0]+[1]*exp(-[2]*x)",
308  "[0]+[1]*exp(-[2]*x)",
309  "pol1",
310  "[0]+[1]*exp(-[2]*x)"
311  };
312 
313  ptnrbins = pTbinVals_.size()-1;
314  double *ptbins = new double[pTbinVals_.size()];
315  for(unsigned int b=0; b<pTbinVals_.size(); b++) ptbins[b] = pTbinVals_[b];
316  double *etabins;
317  if(objectType_ != "met"){
318  etanrbins = etabinVals_.size()-1;
319  etabins = new double[etabinVals_.size()];
320  for(unsigned int b=0; b<etabinVals_.size(); b++) etabins[b] = etabinVals_[b];
321  }else{
322  etanrbins = 1;
323  etabins = new double[2];
324  etabins[0] = 0; etabins[1] = 5.;
325  }
326 
327 
328  //define the histograms booked
329  for(Int_t ro=0; ro<8; ro++) {
330  for(Int_t etab=0; etab<etanrbins; etab++) {
331  for(Int_t ptb=0; ptb<ptnrbins; ptb++) {
332  TString obsName = objectType_; obsName += resObsName[ro]; obsName += "_etabin"; obsName += etab; obsName += "_ptbin";
333  obsName += ptb;
334  hResPtEtaBin[ro][etab][ptb] = fs->make<TH1F>(obsName,obsName,resObsNrBins,resObsMin[ro],resObsMax[ro]);
335  fResPtEtaBin[ro][etab][ptb] = fs->make<TF1>("F_"+obsName,"gaus");
336  }
337  TString obsName2 = objectType_; obsName2 += resObsName[ro]; obsName2 += "_etabin"; obsName2 += etab;
338  hResEtaBin[ro][etab] = fs->make<TH1F>(obsName2,obsName2,ptnrbins,ptbins);
339  fResEtaBin[ro][etab] = fs->make<TF1>("F_"+obsName2,resObsVsPtFit[ro],pTbinVals_[0],pTbinVals_[pTbinVals_.size()-1]);
340  }
341  }
342  tResVar = fs->make< TTree >("tResVar","Resolution tree");
343 
344  delete [] etabins;
345  delete [] ptbins;
346 
347 }
348 
349 // ------------ method called once each job just after ending the event loop ------------
350 void
351 ResolutionCreator::endJob() {
352  TString resObsName2[8] = {"a","b","c","d","theta","phi","et","eta"};
353  Int_t ro=0;
354  Double_t pt=0.;
355  Double_t eta=0.;
356  Double_t value,error;
357 
358  tResVar->Branch("Pt",&pt,"Pt/D");
359  tResVar->Branch("Eta",&eta,"Eta/D");
360  tResVar->Branch("ro",&ro,"ro/I");
361  tResVar->Branch("value",&value,"value/D");
362  tResVar->Branch("error",&error,"error/D");
363 
364  for(ro=0; ro<8; ro++) {
365  for(int etab=0; etab<etanrbins; etab++) {
366  //CD set eta at the center of the bin
367  eta = etanrbins > 1 ? (etabinVals_[etab]+etabinVals_[etab+1])/2. : 2.5 ;
368  for(int ptb=0; ptb<ptnrbins; ptb++) {
369  //CD set pt at the center of the bin
370  pt = (pTbinVals_[ptb]+pTbinVals_[ptb+1])/2.;
371  double maxcontent = 0.;
372  int maxbin = 0;
373  for(int nb=1; nb<hResPtEtaBin[ro][etab][ptb]->GetNbinsX(); nb ++){
374  if (hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb)>maxcontent) {
375  maxcontent = hResPtEtaBin[ro][etab][ptb]->GetBinContent(nb);
376  maxbin = nb;
377  }
378  }
379  int range = (int)(hResPtEtaBin[ro][etab][ptb]->GetNbinsX()/6); //in order that ~1/3 of X-axis range is fitted
380  fResPtEtaBin[ro][etab][ptb] -> SetRange(hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin-range),
381  hResPtEtaBin[ro][etab][ptb]->GetBinCenter(maxbin+range));
382  fResPtEtaBin[ro][etab][ptb] -> SetParameters(hResPtEtaBin[ro][etab][ptb] -> GetMaximum(),
383  hResPtEtaBin[ro][etab][ptb] -> GetMean(),
384  hResPtEtaBin[ro][etab][ptb] -> GetRMS());
385  hResPtEtaBin[ro][etab][ptb] -> Fit(fResPtEtaBin[ro][etab][ptb]->GetName(),"RQ");
386  hResEtaBin[ro][etab] -> SetBinContent(ptb+1,fResPtEtaBin[ro][etab][ptb]->GetParameter(2));
387  hResEtaBin[ro][etab] -> SetBinError(ptb+1,fResPtEtaBin[ro][etab][ptb]->GetParError(2));
388  //CD: Fill the tree
389  value = fResPtEtaBin[ro][etab][ptb]->GetParameter(2); //parameter value
390  error = fResPtEtaBin[ro][etab][ptb]->GetParError(2); //parameter error
391  tResVar->Fill();
392  }
393  //CD: add a fake entry in pt=0 for the NN training
394  // for that, use a linear extrapolation.
395  pt = 0.;
396  value = ((pTbinVals_[0]+pTbinVals_[1])/2.)*(fResPtEtaBin[ro][etab][0]->GetParameter(2)-fResPtEtaBin[ro][etab][1]->GetParameter(2))/((pTbinVals_[2]-pTbinVals_[0])/2.)+fResPtEtaBin[ro][etab][0]->GetParameter(2);
397  error = fResPtEtaBin[ro][etab][0]->GetParError(2)+fResPtEtaBin[ro][etab][1]->GetParError(2);
398  tResVar->Fill();
399  // standard fit
400  hResEtaBin[ro][etab] -> Fit(fResEtaBin[ro][etab]->GetName(),"RQ");
401  }
402  }
403  if(objectType_ == "lJets" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for light jets \n";
404  if(objectType_ == "bJets" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for bjets \n";
405  if(objectType_ == "muon" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for muons \n";
406  if(objectType_ == "electron" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for electrons \n";
407  if(objectType_ == "tau" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for taus \n";
408  if(objectType_ == "met" && nrFilled == 0) edm::LogProblem ("SummaryError") << "No plots filled for met \n";
409 
410  edm::LogVerbatim ("MainResults") << " \n\n";
411  edm::LogVerbatim ("MainResults") << " ----------------------------------------------";
412  edm::LogVerbatim ("MainResults") << " ----------------------------------------------";
413  edm::LogVerbatim ("MainResults") << " Resolutions on "<< objectType_ << " with nrfilled: "<<nrFilled;
414  edm::LogVerbatim ("MainResults") << " ----------------------------------------------";
415  edm::LogVerbatim ("MainResults") << " ----------------------------------------------";
416  if(nrFilled != 0 && objectType_ != "met") {
417  for(ro=0; ro<8; ro++) {
418  edm::LogVerbatim ("MainResults") << "-------------------- ";
419  edm::LogVerbatim ("MainResults") << "\n Resolutions on " << resObsName2[ro] << "\n";
420  edm::LogVerbatim ("MainResults") << "-------------------- ";
421  for(int etab=0; etab<etanrbins; etab++) {
422  if(nrFilled != 0 && ro != 6) {
423  if(etab == 0){
424  edm::LogVerbatim ("MainResults") << "if(fabs(eta)<"<<etabinVals_[etab+1] <<") res = " <<
425  fResEtaBin[ro][etab]->GetParameter(0) << "+" << fResEtaBin[ro][etab]->GetParameter(1)
426  << "*exp(-(" <<fResEtaBin[ro][etab]->GetParameter(2) << "*pt));";
427  }else{
428  edm::LogVerbatim ("MainResults") << "else if(fabs(eta)<"<<etabinVals_[etab+1] <<") res = " <<
429  fResEtaBin[ro][etab]->GetParameter(0) << "+" << fResEtaBin[ro][etab]->GetParameter(1)
430  << "*exp(-(" <<fResEtaBin[ro][etab]->GetParameter(2) << "*pt));";
431  }
432  }else if(nrFilled != 0 && ro == 6){
433  if(etab == 0){
434  edm::LogVerbatim ("MainResults") << "if(fabs(eta)<"<<etabinVals_[etab+1] <<") res = " <<
435  fResEtaBin[ro][etab]->GetParameter(0) << "+" << fResEtaBin[ro][etab]->GetParameter(1)
436  << "*pt;";
437  }else{
438  edm::LogVerbatim ("MainResults") << "else if(fabs(eta)<"<<etabinVals_[etab+1] <<") res = " <<
439  fResEtaBin[ro][etab]->GetParameter(0) << "+" << fResEtaBin[ro][etab]->GetParameter(1)
440  << "*pt;";
441 
442  }
443  }
444  }
445  }
446  }else if(nrFilled != 0 && objectType_ == "met"){
447  for(ro=0; ro<8; ro++) {
448  edm::LogVerbatim ("MainResults") << "-------------------- ";
449  edm::LogVerbatim ("MainResults") << "\n Resolutions on " << resObsName2[ro] << "\n";
450  edm::LogVerbatim ("MainResults") << "-------------------- ";
451  if(nrFilled != 0 && ro != 6) {
452  edm::LogVerbatim ("MainResults") << "res = " <<
453  fResEtaBin[ro][0]->GetParameter(0) << "+" << fResEtaBin[ro][0]->GetParameter(1)
454  << "*exp(-(" <<fResEtaBin[ro][0]->GetParameter(2) << "*pt));";
455  }else if(nrFilled != 0 && ro == 6){
456  edm::LogVerbatim ("MainResults") << "res = " <<
457  fResEtaBin[ro][0]->GetParameter(0) << "+" << fResEtaBin[ro][0]->GetParameter(1) << "*pt;";
458  }
459  }
460  }
461 }
462 
463 //define this as a plug-in
464 DEFINE_FWK_MODULE(ResolutionCreator);
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:161
TtGenEvent::isSemiLeptonic
bool isSemiLeptonic(bool excludeTauLeptons=false) const
check if the event can be classified as semi-laptonic
Definition: TtGenEvent.h:38
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:579
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
ResolutionCreator::endJob
void endJob() override
Definition: ResolutionCreator.cc:351
ResolutionCreator::beginJob
void beginJob() override
Definition: ResolutionCreator.cc:276
reco::CompositeRefCandidateT::numberOfMothers
size_t numberOfMothers() const override
number of mothers
ResolutionCreator::fResEtaBin
TF1 * fResEtaBin[10][20]
Definition: ResolutionCreator.cc:67
reco::Candidate::mother
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
hitfit::phidiff
double phidiff(double phi)
Normalized difference in azimuthal angles to a range between .
Definition: fourvec.cc:230
ResolutionCreator::genEvtToken_
edm::EDGetTokenT< TtGenEvent > genEvtToken_
Definition: ResolutionCreator.cc:49
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
bgen
double bgen
Definition: HydjetWrapper.h:47
ResolutionCreator::jetsToken_
edm::EDGetTokenT< std::vector< pat::Jet > > jetsToken_
Definition: ResolutionCreator.cc:53
iEvent
int iEvent
Definition: GenABIO.cc:230
ResolutionCreator::tausToken_
edm::EDGetTokenT< std::vector< pat::Tau > > tausToken_
Definition: ResolutionCreator.cc:55
ResolutionCreator::etabinVals_
std::vector< double > etabinVals_
Definition: ResolutionCreator.cc:56
Fit
Definition: Fit.h:34
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
fwrapper::jets
vector< PseudoJet > jets
Definition: fastjetfortran_madfks.cc:42
reco::Candidate::numberOfMothers
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
HiRegitMuonDetachedTripletStep_cff.DeltaR
DeltaR
Definition: HiRegitMuonDetachedTripletStep_cff.py:22
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:110
ResolutionCreator::~ResolutionCreator
~ResolutionCreator() override
Definition: ResolutionCreator.cc:99
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
extraflags_cff.muons
muons
Definition: extraflags_cff.py:17
reco::LeafCandidate::p4
const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:99
ResolutionCreator::hResEtaBin
TH1F * hResEtaBin[10][20]
Definition: ResolutionCreator.cc:69
ResolutionCreator::fResPtEtaBin
TF1 * fResPtEtaBin[10][20][20]
Definition: ResolutionCreator.cc:66
TtGenEvent::singleNeutrino
const reco::GenParticle * singleNeutrino(bool excludeTauLeptons=false) const
return single neutrino if available; 0 else
Definition: TtGenEvent.cc:150
ResolutionCreator::hResPtEtaBin
TH1F * hResPtEtaBin[10][20][20]
Definition: ResolutionCreator.cc:68
ResolutionCreator::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: ResolutionCreator.cc:110
ResolutionCreator::ResolutionCreator
ResolutionCreator(const edm::ParameterSet &)
Definition: ResolutionCreator.cc:77
ResolutionCreator::electronsToken_
edm::EDGetTokenT< std::vector< pat::Electron > > electronsToken_
Definition: ResolutionCreator.cc:51
b
double b
Definition: hdecay.h:120
TopGenEvent::particles
const reco::GenParticleCollection & particles() const
return particles of decay chain
Definition: TopGenEvent.h:53
ResolutionCreator::metsToken_
edm::EDGetTokenT< std::vector< pat::MET > > metsToken_
Definition: ResolutionCreator.cc:54
reco::LeafCandidate::status
int status() const final
status word
Definition: LeafCandidate.h:165
ResolutionCreator::pTbinVals_
std::vector< double > pTbinVals_
Definition: ResolutionCreator.cc:56
TtSemiLepSignalSelMVAComputer_cfi.mets
mets
met input
Definition: TtSemiLepSignalSelMVAComputer_cfi.py:8
nano_cff.taus
taus
Definition: nano_cff.py:34
ResolutionCreator::muonsToken_
edm::EDGetTokenT< std::vector< pat::Muon > > muonsToken_
Definition: ResolutionCreator.cc:52
reco::CompositeRefCandidateT::mother
const Candidate * mother(size_type=0) const override
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
TtGenEvtProducer_cfi.genEvt
genEvt
Definition: TtGenEvtProducer_cfi.py:7
electrons_cff.electrons
electrons
Definition: electrons_cff.py:128
Generated for CMSSW Reference Manual by doxygen 1.8.11