CMS 3D CMS Logo
ErrorsPropagationAnalyzer.cc
Go to the documentation of this file.
1 #ifndef ERRORSANALYZER_CC
2 #define ERRORSANALYZER_CC
3 
4 #include "ErrorsPropagationAnalyzer.h"
5 
6 ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig) :
7  treeFileName_( iConfig.getParameter<std::string>("InputFileName") ),
8  resolFitType_( iConfig.getParameter<int>("ResolFitType") ),
9  maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
10  outputFileName_( iConfig.getParameter<std::string>("OutputFileName") ),
11  ptBins_( iConfig.getParameter<int>("PtBins") ),
12  ptMin_( iConfig.getParameter<double>("PtMin") ),
13  ptMax_( iConfig.getParameter<double>("PtMax") ),
14  etaBins_( iConfig.getParameter<int>("EtaBins") ),
15  etaMin_( iConfig.getParameter<double>("EtaMin") ),
16  etaMax_( iConfig.getParameter<double>("EtaMax") ),
17  debug_( iConfig.getParameter<bool>("Debug") ),
18  ptMinCut_( iConfig.getUntrackedParameter<double>("PtMinCut", 0.) ),
19  ptMaxCut_( iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.) ),
20  etaMinCut_( iConfig.getUntrackedParameter<double>("EtaMinCut", 0.) ),
21  etaMaxCut_( iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.) )
22 {
23  parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
24  errors_ = iConfig.getParameter<std::vector<double> >("Errors");
25  errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
26 
27  if( (parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size()) ) {
28  std::cout << "Error: parameters and errors have different number of values" << std::endl;
29  exit(1);
30  }
31 
32  fillValueError();
33 
34  sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
35  sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
36  sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
37 
38  sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
39  sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
40  sigmaPtVsEtaMinusErr_ = new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
41 
42  sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
43  sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
44  sigmaMassVsPtMinusErr_ = new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
45 
46  sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
47  sigmaMassVsEtaPlusErr_ = new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
48  sigmaMassVsEtaMinusErr_ = new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
49 
50  sigmaMassOverMassVsPt_ = new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
51  sigmaMassOverMassVsPtPlusErr_ = new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
52  sigmaMassOverMassVsPtMinusErr_ = new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
53 
54  sigmaMassOverMassVsEta_ = new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
55  sigmaMassOverMassVsEtaPlusErr_ = new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
56  sigmaMassOverMassVsEtaMinusErr_ = new TProfile("sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
57 
58  sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
59  sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
60 }
61 
62 void ErrorsPropagationAnalyzer::fillValueError()
63 {
64  valuePlusError_.resize(parameters_.size());
65  valueMinusError_.resize(parameters_.size());
66 
67  std::vector<double>::const_iterator parIt = parameters_.begin();
68  std::vector<double>::const_iterator errIt = errors_.begin();
69  std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
70  int i=0;
71  for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
72  valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
73  valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
74  }
75 }
76 
77 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer()
78 {
79  gROOT->SetStyle("Plain");
80 
81  fillHistograms();
82 
83  TFile * outputFile = new TFile(outputFileName_, "RECREATE");
84 
85  drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
86  drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");
87 
88  drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
89  drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");
90 
91  drawHistograms(sigmaMassOverMassVsEta_, sigmaMassOverMassVsEtaPlusErr_, sigmaMassOverMassVsEtaMinusErr_, "sigmaMassOverMassVsEta", "#sigma(M)/M");
92  drawHistograms(sigmaMassOverMassVsPt_, sigmaMassOverMassVsPtPlusErr_, sigmaMassOverMassVsPtMinusErr_, "sigmaMassOverMassVsPt", "#sigma(M)/M");
93 
94  sigmaPtVsPtDiff_->Write();
95  sigmaPtVsEtaDiff_->Write();
96 
97  outputFile->Write();
98  outputFile->Close();
99 }
100 
101 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile * histo, const TProfile * histoPlusErr,
102  const TProfile * histoMinusErr, const TString & type,
103  const TString & yLabel)
104 {
105  TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
106  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
107  TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
108  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
109  TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
110  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
111 
112  TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
113  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
114  TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
115  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
116  TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
117  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
118 
119  TH1D * sigmaMassOverMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
120  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
121  TH1D * sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
122  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
123  TH1D * sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
124  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
125 
126  TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
127  for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
128  sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
129  sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
130  sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
131  }
132  int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
133  // Draw TGraph with asymmetric errors
134  Double_t * values = sigmaPtVsEtaTH1D->GetArray();
135  Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
136  Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
137  double * posErrors = new double[numBins];
138  double * negErrors = new double[numBins];
139 
140  TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
141  TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);
142 
143  for( int i=1; i<=numBins; ++i ) {
144  // std::cout << "filling " << i << std::endl;
145  posErrors[i-1] = valuesPlus[i] - values[i];
146  if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
147  else negErrors[i-1] = values[i] - valuesMinus[i];
148 
149  graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
150  graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
151  }
152 
153  canvas->Draw();
154 
155  graphAsymmErrors->SetTitle("");
156  graphAsymmErrors->SetFillColor(kGray);
157  graphAsymmErrors->Draw("A2");
158  TString title(type);
159  if( type == "Eta" ) title = "#eta";
160  graphAsymmErrors->GetXaxis()->SetTitle(title);
161  graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
162  graph->Draw();
163 
164  // graph->SetLineColor(kGray);
165  // graph->SetMarkerColor(kBlack);
166  // graph->Draw("AP");
167 
168 // if( debug_ ) {
169 // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
170 // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
171 // }
172 // else {
173 // sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
174 // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
175 // sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
176 // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
177 // }
178 // sigmaPtVsEtaPlusErrTH1D->Draw();
179 // sigmaPtVsEtaTH1D->Draw("SAME");
180 // sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
181 
182  sigmaPtVsEtaPlusErrTH1D->Write();
183  sigmaPtVsEtaTH1D->Write();
184  sigmaPtVsEtaMinusErrTH1D->Write();
185 
186  sigmaPtVsEtaPlusErr_->Write();
187  sigmaPtVsEta_->Write();
188  sigmaPtVsEtaMinusErr_->Write();
189 
190  // Mass
191  sigmaMassVsEtaPlusErrTH1D->Write();
192  sigmaMassVsEtaTH1D->Write();
193  sigmaMassVsEtaMinusErrTH1D->Write();
194 
195  sigmaMassOverMassVsEtaPlusErrTH1D->Write();
196  sigmaMassOverMassVsEtaTH1D->Write();
197  sigmaMassOverMassVsEtaMinusErrTH1D->Write();
198 
199  sigmaMassVsEtaPlusErr_->Write();
200  sigmaMassVsEta_->Write();
201  sigmaMassVsEtaMinusErr_->Write();
202 
203  sigmaMassOverMassVsEtaPlusErr_->Write();
204  sigmaMassOverMassVsEta_->Write();
205  sigmaMassOverMassVsEtaMinusErr_->Write();
206 
207  canvas->Write();
208 }
209 
210 // ------------ method called to for each event ------------
211 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
212 {
213 }
214 
215 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
216  const lorentzVector& mu2,
217  const std::vector< double >& parval,
218  const double& sigmaPt1,
219  const double& sigmaPt2)
220 {
221  double * p = new double[(int)(parval.size())];
222  std::vector<double>::const_iterator it = parval.begin();
223  int id = 0;
224  for ( ; it!=parval.end(); ++it, ++id) {
225  p[id] = *it;
226  }
227  double massRes = massResolution (mu1, mu2, p, sigmaPt1, sigmaPt2);
228  delete[] p;
229  return massRes;
230 }
231 
232 double ErrorsPropagationAnalyzer::massResolution( const lorentzVector& mu1,
233  const lorentzVector& mu2,
234  double* parval,
235  const double & sigmaPt1,
236  const double & sigmaPt2)
237 {
238  double mass = (mu1+mu2).mass();
239  double pt1 = mu1.Pt();
240  double phi1 = mu1.Phi();
241  double eta1 = mu1.Eta();
242  double theta1 = 2*atan(exp(-eta1));
243  double pt2 = mu2.Pt();
244  double phi2 = mu2.Phi();
245  double eta2 = mu2.Eta();
246  double theta2 = 2*atan(exp(-eta2));
247 
248  // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
249  // -----------------------------------------------------------------------------------------------------------
250  double dmdpt1 = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)/(std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2))-
251  pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
252  double dmdpt2 = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)/(std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2))-
253  pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
254  double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
255  double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
256  double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
257  sqrt((std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)/(std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)) -
258  pt1*pt2*cos(theta2)/sin(theta2))/mass;
259  double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
260  sqrt((std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)/(std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)) -
261  pt2*pt1*cos(theta1)/sin(theta1))/mass;
262 
263  // Resolution parameters:
264  // ----------------------
265  // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
266  // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
267  double sigma_pt1 = sigmaPt1;
268  double sigma_pt2 = sigmaPt2;
269  double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt1,eta1,parval );
270  double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt2,eta2,parval );
271  double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt1,eta1,parval );
272  double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt2,eta2,parval );
273  double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2( pt1, eta1, pt2, eta2, parval );
274 
275  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
276  // multiply it by pt.
277  double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2)+
278  std::pow(dmdphi1*sigma_phi1,2)+std::pow(dmdphi2*sigma_phi2,2)+
279  std::pow(dmdcotgth1*sigma_cotgth1,2)+std::pow(dmdcotgth2*sigma_cotgth2,2)+
280  2*dmdpt1*dmdpt2*cov_pt1pt2);
281 
282  return mass_res;
283 }
284 
285 
286 void ErrorsPropagationAnalyzer::fillHistograms()
287 {
288  std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
289  RootTreeHandler rootTreeHandler;
290 
291  typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
292  MuonPairVector savedPair;
293  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
294  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
295  // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
296 
297  resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec = resolutionFunctionVecService( resolFitType_ );
298  MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType_ );
299  MuScleFitUtils::debugMassResol_ = false;
300  MuScleFitUtils::debug = 0;
301  MuScleFitUtils::resfind = std::vector<int>(6, 0);
302 
303  SigmaPt sigmaPt( parameters_, errors_ );
304  SigmaPtDiff sigmaPtDiff;
305 
306  // Loop on all the pairs
307  unsigned int i = 0;
308  MuonPairVector::iterator it = savedPair.begin();
309  std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
310  for( ; it != savedPair.end(); ++it, ++i ) {
311  double pt1 = it->first.pt();
312  double eta1 = it->first.eta();
313  double pt2 = it->second.pt();
314  double eta2 = it->second.eta();
315 
316  if( debug_ ) {
317  std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
318  }
319  // double fabsEta1 = fabs(eta1);
320  // double fabsEta2 = fabs(eta2);
321 
322  if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;
323 
324 
325  // double sigmaPt1 = sigmaPt( eta1 );
326  // double sigmaPt2 = sigmaPt( eta2 );
327 
328  // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
329  // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
330  // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
331  // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
332 
333  double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
334  double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
335  double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_);
336  double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
337  double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_ );
338  double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
339 
340  // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
341  // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
342  // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
343  double sigmaMass = massResolution( it->first, it->second, parameters_, sigmaPt1, sigmaPt2 );
344  double sigmaMassPlusErr = massResolution( it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2 );
345  double sigmaMassMinusErr = massResolution( it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2 );
346 
347  double mass = (it->first + it->second).mass();
348 
349  if( debug_ ) {
350  std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
351  std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
352  std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
353  }
354 
355  // Protections from nans
356  if( pt1 != pt1 ) continue;
357  if( pt2 != pt2 ) continue;
358  if( eta1 != eta1 ) continue;
359  if( eta2 != eta2 ) continue;
360  if( sigmaPt1 != sigmaPt1 ) continue;
361  if( sigmaPt2 != sigmaPt2 ) continue;
362  if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
363  if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
364  if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
365  if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
366  if( sigmaMass != sigmaMass ) continue;
367  if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
368  if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
369  if( mass == 0 ) continue;
370 
371  std::cout << "Muon pair number " << i << std::endl;
372 
373  // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
374  // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
375  // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
376 
377 
378  if( (pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_) ) {
379  sigmaPtVsPt_->Fill(pt1, sigmaPt1);
380  sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
381  sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
382  sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
383  sigmaMassVsPt_->Fill(pt1, sigmaMass*mass);
384  sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr*mass);
385  sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr*mass);
386  sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
387  sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
388  sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
389 
390  sigmaPtVsEta_->Fill(eta1, sigmaPt1);
391  sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
392  sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
393  sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
394  sigmaMassVsEta_->Fill(eta1, sigmaMass*mass);
395  sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr*mass);
396  sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr*mass);
397  sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
398  sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
399  sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
400  }
401  if( (pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_) ) {
402  sigmaPtVsPt_->Fill(pt2, sigmaPt2);
403  sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
404  sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
405  sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
406  sigmaMassVsPt_->Fill(pt2, sigmaMass*mass);
407  sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr*mass);
408  sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr*mass);
409  sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
410  sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
411  sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
412 
413  sigmaPtVsEta_->Fill(eta2, sigmaPt2);
414  sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
415  sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
416  sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
417  sigmaMassVsEta_->Fill(eta2, sigmaMass*mass);
418  sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr*mass);
419  sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr*mass);
420  sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
421  sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
422  sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
423  }
424  }
425 }
426 
427 //define this as a plug-in
428 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
429 
430 #endif
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
mps_fire.i
i
Definition: mps_fire.py:156
SigmaPtDiff
Returns ( sigmaPt/Pt(data)^2 - sigmaPt/Pt(MC)^2 )
Definition: SigmaPtDiff.h:71
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
ErrorsPropagationAnalyzer::valueMinusError_
std::vector< double > valueMinusError_
Definition: ErrorsPropagationAnalyzer.h:93
MuScleFitUtils::debugMassResol_
static bool debugMassResol_
Definition: MuScleFitUtils.h:251
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MuScleFitUtils::debug
static int debug
Definition: MuScleFitUtils.h:113
ErrorsPropagationAnalyzer::valuePlusError_
std::vector< double > valuePlusError_
Definition: ErrorsPropagationAnalyzer.h:92
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
resolutionFunctionBase::sigmaPt
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
iEvent
int iEvent
Definition: GenABIO.cc:230
resolutionFunctionVecService
resolutionFunctionBase< std::vector< double > > * resolutionFunctionVecService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier when receiving a std::...
Definition: Functions.cc:50
hiDetachedQuadStep_cff.pt2
pt2
Definition: hiDetachedQuadStep_cff.py:70
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
resolutionFunctionBase::covPt1Pt2
virtual double covPt1Pt2(const double &pt1, const double &eta1, const double &pt2, const double &eta2, const T &parval)
Definition: Functions.h:685
MuonPairVector
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
Definition: RootTreeHandler.h:13
hiDetachedQuadStep_cff.pt1
pt1
Definition: hiDetachedQuadStep_cff.py:70
SigmaPtDiff::squaredDiff
double squaredDiff(const double &eta)
Definition: SigmaPtDiff.h:139
ErrorsPropagationAnalyzer::massResolution
double massResolution(const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval, const double &sigmaPt1, const double &sigmaPt2)
Modified method to take into account the error.
Definition: ErrorsPropagationAnalyzer.cc:215
resolutionFunctionBase::sigmaCotgTh
virtual double sigmaCotgTh(const double &pt, const double &eta, const T &parval)=0
download_sqlite_cfg.outputFile
outputFile
Definition: download_sqlite_cfg.py:5
ErrorsPropagationAnalyzer::drawHistograms
void drawHistograms(const TProfile *histo, const TProfile *histoPlusErr, const TProfile *histoMinusErr, const TString &type, const TString &yLabel)
Definition: ErrorsPropagationAnalyzer.cc:101
MuScleFitUtils::mMu2
static const double mMu2
Definition: MuScleFitUtils.h:122
fftjetcommon_cfi.title
title
Definition: fftjetcommon_cfi.py:32
MuonErrorMatrixValues_cff.values
values
Definition: MuonErrorMatrixValues_cff.py:6
resolutionFunctionService
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:38
svgfig.canvas
def canvas(sub, attr)
Definition: svgfig.py:481
gather_cfg.cout
cout
Definition: gather_cfg.py:145
ErrorsPropagationAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: ErrorsPropagationAnalyzer.cc:211
ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer
ErrorsPropagationAnalyzer(const edm::ParameterSet &)
Definition: ErrorsPropagationAnalyzer.cc:6
resolutionFunctionBase::sigmaPhi
virtual double sigmaPhi(const double &pt, const double &eta, const T &parval)=0
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
MuScleFitUtils::resfind
static std::vector< int > resfind
Definition: MuScleFitUtils.h:193
resolutionFunctionBase::sigmaPtError
virtual double sigmaPtError(const double &pt, const double &eta, const T &parval, const T &parError)
Definition: Functions.h:679
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
hcalTTPDigis_cfi.id
id
Definition: hcalTTPDigis_cfi.py:7
MuScleFitUtils::resolutionFunction
static resolutionFunctionBase< double * > * resolutionFunction
Definition: MuScleFitUtils.h:136
RootTreeHandler::readTree
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, std::vector< std::pair< unsigned int, unsigned long long > > *evtRun, MuonPairVector *genPair=0)
Definition: RootTreeHandler.h:88
Generated for CMSSW Reference Manual by doxygen 1.8.11