Reference Manual

Go to the documentation of this file.
00001 #ifndef ERRORSANALYZER_CC
00002 #define ERRORSANALYZER_CC
00003 
00004 #include "ErrorsPropagationAnalyzer.h"
00005 
00006 ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig) :
00007   treeFileName_( iConfig.getParameter<std::string>("InputFileName") ),
00008   resolFitType_( iConfig.getParameter<int>("ResolFitType") ),
00009   maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
00010   outputFileName_( iConfig.getParameter<std::string>("OutputFileName") ),
00011   ptBins_( iConfig.getParameter<int>("PtBins") ),
00012   ptMin_( iConfig.getParameter<double>("PtMin") ),
00013   ptMax_( iConfig.getParameter<double>("PtMax") ),
00014   etaBins_( iConfig.getParameter<int>("EtaBins") ),
00015   etaMin_( iConfig.getParameter<double>("EtaMin") ),
00016   etaMax_( iConfig.getParameter<double>("EtaMax") ),
00017   debug_( iConfig.getParameter<bool>("Debug") ),
00018   ptMinCut_( iConfig.getUntrackedParameter<double>("PtMinCut", 0.) ),
00019   ptMaxCut_( iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.) ),
00020   etaMinCut_( iConfig.getUntrackedParameter<double>("EtaMinCut", 0.) ),
00021   etaMaxCut_( iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.) )
00022 {
00023   parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
00024   errors_ = iConfig.getParameter<std::vector<double> >("Errors");
00025   errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
00026 
00027   if( (parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size()) ) {
00028     std::cout << "Error: parameters and errors have different number of values" << std::endl;
00029     exit(1);
00030   }
00031 
00032   fillValueError();
00033 
00034   sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
00035   sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
00036   sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
00037 
00038   sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
00039   sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
00040   sigmaPtVsEtaMinusErr_ = new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
00041 
00042   sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
00043   sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
00044   sigmaMassVsPtMinusErr_ = new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
00045 
00046   sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
00047   sigmaMassVsEtaPlusErr_ = new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
00048   sigmaMassVsEtaMinusErr_ = new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
00049 
00050   sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
00051   sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
00052 }
00053 
00054 void ErrorsPropagationAnalyzer::fillValueError()
00055 {
00056   valuePlusError_.resize(parameters_.size());
00057   valueMinusError_.resize(parameters_.size());
00058 
00059   std::vector<double>::const_iterator parIt = parameters_.begin();
00060   std::vector<double>::const_iterator errIt = errors_.begin();
00061   std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
00062   int i=0;
00063   for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
00064     valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
00065     valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
00066   }
00067 }
00068 
00069 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer()
00070 {
00071   gROOT->SetStyle("Plain");
00072 
00073   fillHistograms();
00074 
00075   TFile * outputFile = new TFile(outputFileName_, "RECREATE");
00076 
00077   drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta");
00078   drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt");
00079 
00080   drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta");
00081   drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt");
00082 
00083   sigmaPtVsPtDiff_->Write();
00084   sigmaPtVsEtaDiff_->Write();  
00085 
00086   outputFile->Write();
00087   outputFile->Close();
00088 }
00089 
00090 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile * histo, const TProfile * histoPlusErr,
00091                                     const TProfile * histoMinusErr, const TString & type)
00092 {
00093   TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
00094                                      histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00095   TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
00096                                             histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00097   TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
00098                                              histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00099 
00100   TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
00101                                        histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00102   TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
00103                                               histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00104   TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
00105                                                histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00106 
00107   TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
00108   for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
00109     sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
00110     sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
00111     sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
00112   }
00113   int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
00114   // Draw TGraph with asymmetric errors
00115   Double_t * values = sigmaPtVsEtaTH1D->GetArray();
00116   Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
00117   Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
00118   double * posErrors = new double[numBins];
00119   double * negErrors = new double[numBins];
00120 
00121   TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
00122   TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);
00123 
00124   for( int i=1; i<=numBins; ++i ) {
00125     // std::cout << "filling " << i << std::endl;
00126     posErrors[i-1] = valuesPlus[i] - values[i];
00127     if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
00128     else negErrors[i-1] = values[i] - valuesMinus[i];
00129 
00130     graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
00131     graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
00132   }
00133 
00134   canvas->Draw();
00135 
00136   graphAsymmErrors->SetTitle("");
00137   graphAsymmErrors->SetFillColor(kGray);
00138   graphAsymmErrors->Draw("A2");
00139   TString title(type);
00140   if( type == "Eta" ) title = "#eta";
00141   graphAsymmErrors->GetXaxis()->SetTitle(title);
00142   graphAsymmErrors->GetYaxis()->SetTitle("#sigmaPt/Pt");
00143   graph->Draw();
00144 
00145   //  graph->SetLineColor(kGray);
00146   //  graph->SetMarkerColor(kBlack);
00147   //  graph->Draw("AP");
00148 
00149 //   if( debug_ ) {
00150 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
00151 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
00152 //   }
00153 //   else {
00154 //     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
00155 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
00156 //     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
00157 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
00158 //   }
00159 //   sigmaPtVsEtaPlusErrTH1D->Draw();
00160 //   sigmaPtVsEtaTH1D->Draw("SAME");
00161 //   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
00162 
00163   sigmaPtVsEtaPlusErrTH1D->Write();
00164   sigmaPtVsEtaTH1D->Write();
00165   sigmaPtVsEtaMinusErrTH1D->Write();
00166 
00167   sigmaPtVsEtaPlusErr_->Write();
00168   sigmaPtVsEta_->Write();
00169   sigmaPtVsEtaMinusErr_->Write();
00170 
00171   // Mass
00172   sigmaMassVsEtaPlusErrTH1D->Write();
00173   sigmaMassVsEtaTH1D->Write();
00174   sigmaMassVsEtaMinusErrTH1D->Write();
00175 
00176   sigmaMassVsEtaPlusErr_->Write();
00177   sigmaMassVsEta_->Write();
00178   sigmaMassVsEtaMinusErr_->Write();
00179 
00180   canvas->Write();
00181 }
00182 
00183 // ------------ method called to for each event  ------------
00184 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00185 {
00186 }
00187 
00188 void ErrorsPropagationAnalyzer::fillHistograms()
00189 {
00190   std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
00191   RootTreeHandler rootTreeHandler;
00192 
00193   typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
00194   MuonPairVector savedPair;
00195   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0);
00196   // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
00197 
00198   resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec = resolutionFunctionVecService( resolFitType_ );
00199   MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType_ );
00200   MuScleFitUtils::debugMassResol_ = false;
00201   MuScleFitUtils::debug = 0;
00202   MuScleFitUtils::resfind = std::vector<int>(6, 0);
00203 
00204   SigmaPt sigmaPt( parameters_, errors_ );
00205   SigmaPtDiff sigmaPtDiff;
00206 
00207   // Loop on all the pairs
00208   unsigned int i = 0;
00209   MuonPairVector::iterator it = savedPair.begin();
00210   std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
00211   for( ; it != savedPair.end(); ++it, ++i ) {
00212     double pt1 = it->first.pt();
00213     double eta1 = it->first.eta();
00214     double pt2 = it->second.pt();
00215     double eta2 = it->second.eta();
00216 
00217     if( debug_ ) {
00218       std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
00219     }
00220     // double fabsEta1 = fabs(eta1);
00221     // double fabsEta2 = fabs(eta2);
00222 
00223     if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;
00224 
00225 
00226     // double sigmaPt1 = sigmaPt( eta1 );
00227     // double sigmaPt2 = sigmaPt( eta2 );
00228 
00229     // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
00230     // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
00231     // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
00232     // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
00233 
00234     double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
00235     double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
00236     double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_);
00237     double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
00238     double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_ );
00239     double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
00240 
00241     double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
00242     double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
00243     double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
00244 
00245     if( debug_ ) {
00246       std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
00247       std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
00248       std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
00249     }
00250 
00251     // Protections from nans
00252     if( pt1 != pt1 ) continue;
00253     if( pt2 != pt2 ) continue;
00254     if( eta1 != eta1 ) continue;
00255     if( eta2 != eta2 ) continue;
00256     if( sigmaPt1 != sigmaPt1 ) continue;
00257     if( sigmaPt2 != sigmaPt2 ) continue;
00258     if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
00259     if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
00260     if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
00261     if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
00262     if( sigmaMass != sigmaMass ) continue;
00263     if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
00264     if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
00265 
00266     std::cout << "Muon pair number " << i << std::endl;
00267 
00268     // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
00269     // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
00270     // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
00271 
00272 
00273     if( (pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_) ) {
00274       sigmaPtVsPt_->Fill(pt1, sigmaPt1);
00275       sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
00276       sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
00277       sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
00278       sigmaMassVsPt_->Fill(pt1, sigmaMass);
00279       sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
00280       sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
00281 
00282       sigmaPtVsEta_->Fill(eta1, sigmaPt1);
00283       sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
00284       sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
00285       sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
00286       sigmaMassVsEta_->Fill(eta1, sigmaMass);
00287       sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
00288       sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
00289     }
00290     if( (pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_) ) {
00291       sigmaPtVsPt_->Fill(pt2, sigmaPt2);
00292       sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
00293       sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
00294       sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
00295       sigmaMassVsPt_->Fill(pt2, sigmaMass);
00296       sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
00297       sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
00298 
00299       sigmaPtVsEta_->Fill(eta2, sigmaPt2);
00300       sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
00301       sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
00302       sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
00303       sigmaMassVsEta_->Fill(eta2, sigmaMass);
00304       sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
00305       sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
00306     }
00307   }
00308 }
00309 
00310 //define this as a plug-in
00311 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
00312 
00313 #endif
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/MuonAnalysis/MomentumScaleCalibration/plugins/ErrorsPropagationAnalyzer.cc
