Reference Manual

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