CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ErrorsPropagationAnalyzer.cc
Go to the documentation of this file.
1 #ifndef ERRORSANALYZER_CC
2 #define ERRORSANALYZER_CC
3 
5 
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 
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  sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
51  sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
52 }
53 
55 {
56  valuePlusError_.resize(parameters_.size());
57  valueMinusError_.resize(parameters_.size());
58 
59  std::vector<double>::const_iterator parIt = parameters_.begin();
60  std::vector<double>::const_iterator errIt = errors_.begin();
61  std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
62  int i=0;
63  for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
64  valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
65  valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
66  }
67 }
68 
70 {
71  gROOT->SetStyle("Plain");
72 
74 
75  TFile * outputFile = new TFile(outputFileName_, "RECREATE");
76 
79 
82 
83  sigmaPtVsPtDiff_->Write();
84  sigmaPtVsEtaDiff_->Write();
85 
86  outputFile->Write();
87  outputFile->Close();
88 }
89 
90 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile * histo, const TProfile * histoPlusErr,
91  const TProfile * histoMinusErr, const TString & type)
92 {
93  TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
94  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
95  TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
96  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
97  TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
98  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
99 
100  TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
101  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
102  TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
103  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
104  TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
105  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
106 
107  TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
108  for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
109  sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
110  sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
111  sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
112  }
113  int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
114  // Draw TGraph with asymmetric errors
115  Double_t * values = sigmaPtVsEtaTH1D->GetArray();
116  Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
117  Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
118  double * posErrors = new double[numBins];
119  double * negErrors = new double[numBins];
120 
121  TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
122  TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);
123 
124  for( int i=1; i<=numBins; ++i ) {
125  // std::cout << "filling " << i << std::endl;
126  posErrors[i-1] = valuesPlus[i] - values[i];
127  if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
128  else negErrors[i-1] = values[i] - valuesMinus[i];
129 
130  graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
131  graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
132  }
133 
134  canvas->Draw();
135 
136  graphAsymmErrors->SetTitle("");
137  graphAsymmErrors->SetFillColor(kGray);
138  graphAsymmErrors->Draw("A2");
139  TString title(type);
140  if( type == "Eta" ) title = "#eta";
141  graphAsymmErrors->GetXaxis()->SetTitle(title);
142  graphAsymmErrors->GetYaxis()->SetTitle("#sigmaPt/Pt");
143  graph->Draw();
144 
145  // graph->SetLineColor(kGray);
146  // graph->SetMarkerColor(kBlack);
147  // graph->Draw("AP");
148 
149 // if( debug_ ) {
150 // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
151 // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
152 // }
153 // else {
154 // sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
155 // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
156 // sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
157 // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
158 // }
159 // sigmaPtVsEtaPlusErrTH1D->Draw();
160 // sigmaPtVsEtaTH1D->Draw("SAME");
161 // sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
162 
163  sigmaPtVsEtaPlusErrTH1D->Write();
164  sigmaPtVsEtaTH1D->Write();
165  sigmaPtVsEtaMinusErrTH1D->Write();
166 
167  sigmaPtVsEtaPlusErr_->Write();
168  sigmaPtVsEta_->Write();
169  sigmaPtVsEtaMinusErr_->Write();
170 
171  // Mass
172  sigmaMassVsEtaPlusErrTH1D->Write();
173  sigmaMassVsEtaTH1D->Write();
174  sigmaMassVsEtaMinusErrTH1D->Write();
175 
176  sigmaMassVsEtaPlusErr_->Write();
177  sigmaMassVsEta_->Write();
178  sigmaMassVsEtaMinusErr_->Write();
179 
180  canvas->Write();
181 }
182 
183 // ------------ method called to for each event ------------
185 {
186 }
187 
189 {
190  std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
191  RootTreeHandler rootTreeHandler;
192 
193  typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
194  MuonPairVector savedPair;
195  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0);
196  // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
197 
202  MuScleFitUtils::resfind = std::vector<int>(6, 0);
203 
204  SigmaPt sigmaPt( parameters_, errors_ );
205  SigmaPtDiff sigmaPtDiff;
206 
207  // Loop on all the pairs
208  unsigned int i = 0;
209  MuonPairVector::iterator it = savedPair.begin();
210  std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
211  for( ; it != savedPair.end(); ++it, ++i ) {
212  double pt1 = it->first.pt();
213  double eta1 = it->first.eta();
214  double pt2 = it->second.pt();
215  double eta2 = it->second.eta();
216 
217  if( debug_ ) {
218  std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
219  }
220  // double fabsEta1 = fabs(eta1);
221  // double fabsEta2 = fabs(eta2);
222 
223  if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;
224 
225 
226  // double sigmaPt1 = sigmaPt( eta1 );
227  // double sigmaPt2 = sigmaPt( eta2 );
228 
229  // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
230  // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
231  // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
232  // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
233 
234  double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
235  double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
236  double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_);
237  double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
238  double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_ );
239  double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
240 
241  double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
242  double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
243  double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
244 
245  if( debug_ ) {
246  std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
247  std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
248  std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
249  }
250 
251  // Protections from nans
252  if( pt1 != pt1 ) continue;
253  if( pt2 != pt2 ) continue;
254  if( eta1 != eta1 ) continue;
255  if( eta2 != eta2 ) continue;
256  if( sigmaPt1 != sigmaPt1 ) continue;
257  if( sigmaPt2 != sigmaPt2 ) continue;
258  if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
259  if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
260  if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
261  if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
262  if( sigmaMass != sigmaMass ) continue;
263  if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
264  if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
265 
266  std::cout << "Muon pair number " << i << std::endl;
267 
268  // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
269  // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
270  // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
271 
272 
273  if( (pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_) ) {
274  sigmaPtVsPt_->Fill(pt1, sigmaPt1);
275  sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
276  sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
277  sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
278  sigmaMassVsPt_->Fill(pt1, sigmaMass);
279  sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
280  sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
281 
282  sigmaPtVsEta_->Fill(eta1, sigmaPt1);
283  sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
284  sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
285  sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
286  sigmaMassVsEta_->Fill(eta1, sigmaMass);
287  sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
288  sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
289  }
290  if( (pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_) ) {
291  sigmaPtVsPt_->Fill(pt2, sigmaPt2);
292  sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
293  sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
294  sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
295  sigmaMassVsPt_->Fill(pt2, sigmaMass);
296  sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
297  sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
298 
299  sigmaPtVsEta_->Fill(eta2, sigmaPt2);
300  sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
301  sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
302  sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
303  sigmaMassVsEta_->Fill(eta2, sigmaMass);
304  sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
305  sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
306  }
307  }
308 }
309 
310 //define this as a plug-in
312 
313 #endif
type
Definition: HCALResponse.h:22
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
Returns ( sigmaPt/Pt(data)^2 - sigmaPt/Pt(MC)^2 )
Definition: SigmaPtDiff.h:68
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< double > valueMinusError_
static bool debugMassResol_
static int debug
tuple histo
Definition: trackerHits.py:12
std::vector< double > valuePlusError_
void drawHistograms(const TProfile *histo, const TProfile *histoPlusErr, const TProfile *histoMinusErr, const TString &type)
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
int iEvent
Definition: GenABIO.cc:243
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:116
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, MuonPairVector *genPair=0)
double squaredDiff(const double &eta)
Definition: SigmaPtDiff.h:136
Definition: adjgraph.h:15
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:88
tuple cout
Definition: gather_cfg.py:41
virtual void analyze(const edm::Event &, const edm::EventSetup &)
ErrorsPropagationAnalyzer(const edm::ParameterSet &)
static std::vector< int > resfind
virtual double sigmaPtError(const double &pt, const double &eta, const T &parval, const T &parError)
Definition: Functions.h:1525
static resolutionFunctionBase< double * > * resolutionFunction