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   sigmaMassOverMassVsPt_ = new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
00051   sigmaMassOverMassVsPtPlusErr_ = new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
00052   sigmaMassOverMassVsPtMinusErr_ = new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
00053 
00054   sigmaMassOverMassVsEta_ = new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
00055   sigmaMassOverMassVsEtaPlusErr_ = new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
00056   sigmaMassOverMassVsEtaMinusErr_ = new TProfile("sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
00057 
00058   sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
00059   sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
00060 }
00061 
00062 void ErrorsPropagationAnalyzer::fillValueError()
00063 {
00064   valuePlusError_.resize(parameters_.size());
00065   valueMinusError_.resize(parameters_.size());
00066 
00067   std::vector<double>::const_iterator parIt = parameters_.begin();
00068   std::vector<double>::const_iterator errIt = errors_.begin();
00069   std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
00070   int i=0;
00071   for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
00072     valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
00073     valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
00074   }
00075 }
00076 
00077 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer()
00078 {
00079   gROOT->SetStyle("Plain");
00080 
00081   fillHistograms();
00082 
00083   TFile * outputFile = new TFile(outputFileName_, "RECREATE");
00084 
00085   drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
00086   drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");
00087 
00088   drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
00089   drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");
00090 
00091   drawHistograms(sigmaMassOverMassVsEta_, sigmaMassOverMassVsEtaPlusErr_, sigmaMassOverMassVsEtaMinusErr_, "sigmaMassOverMassVsEta", "#sigma(M)/M");
00092   drawHistograms(sigmaMassOverMassVsPt_, sigmaMassOverMassVsPtPlusErr_, sigmaMassOverMassVsPtMinusErr_, "sigmaMassOverMassVsPt", "#sigma(M)/M");
00093 
00094   sigmaPtVsPtDiff_->Write();
00095   sigmaPtVsEtaDiff_->Write();  
00096 
00097   outputFile->Write();
00098   outputFile->Close();
00099 }
00100 
00101 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile * histo, const TProfile * histoPlusErr,
00102                                                const TProfile * histoMinusErr, const TString & type,
00103                                                const TString & yLabel)
00104 {
00105   TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
00106                                      histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00107   TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
00108                                             histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00109   TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
00110                                              histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00111 
00112   TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
00113                                        histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00114   TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
00115                                               histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00116   TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
00117                                                histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00118 
00119   TH1D * sigmaMassOverMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
00120                                                histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00121   TH1D * sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
00122                                                       histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00123   TH1D * sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
00124                                                        histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
00125 
00126   TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
00127   for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
00128     sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
00129     sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
00130     sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
00131   }
00132   int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
00133   // Draw TGraph with asymmetric errors
00134   Double_t * values = sigmaPtVsEtaTH1D->GetArray();
00135   Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
00136   Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
00137   double * posErrors = new double[numBins];
00138   double * negErrors = new double[numBins];
00139 
00140   TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
00141   TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);
00142 
00143   for( int i=1; i<=numBins; ++i ) {
00144     // std::cout << "filling " << i << std::endl;
00145     posErrors[i-1] = valuesPlus[i] - values[i];
00146     if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
00147     else negErrors[i-1] = values[i] - valuesMinus[i];
00148 
00149     graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
00150     graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
00151   }
00152 
00153   canvas->Draw();
00154 
00155   graphAsymmErrors->SetTitle("");
00156   graphAsymmErrors->SetFillColor(kGray);
00157   graphAsymmErrors->Draw("A2");
00158   TString title(type);
00159   if( type == "Eta" ) title = "#eta";
00160   graphAsymmErrors->GetXaxis()->SetTitle(title);
00161   graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
00162   graph->Draw();
00163 
00164   //  graph->SetLineColor(kGray);
00165   //  graph->SetMarkerColor(kBlack);
00166   //  graph->Draw("AP");
00167 
00168 //   if( debug_ ) {
00169 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
00170 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
00171 //   }
00172 //   else {
00173 //     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
00174 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
00175 //     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
00176 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
00177 //   }
00178 //   sigmaPtVsEtaPlusErrTH1D->Draw();
00179 //   sigmaPtVsEtaTH1D->Draw("SAME");
00180 //   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
00181 
00182   sigmaPtVsEtaPlusErrTH1D->Write();
00183   sigmaPtVsEtaTH1D->Write();
00184   sigmaPtVsEtaMinusErrTH1D->Write();
00185 
00186   sigmaPtVsEtaPlusErr_->Write();
00187   sigmaPtVsEta_->Write();
00188   sigmaPtVsEtaMinusErr_->Write();
00189 
00190   // Mass
00191   sigmaMassVsEtaPlusErrTH1D->Write();
00192   sigmaMassVsEtaTH1D->Write();
00193   sigmaMassVsEtaMinusErrTH1D->Write();
00194 
00195   sigmaMassOverMassVsEtaPlusErrTH1D->Write();
00196   sigmaMassOverMassVsEtaTH1D->Write();
00197   sigmaMassOverMassVsEtaMinusErrTH1D->Write();
00198 
00199   sigmaMassVsEtaPlusErr_->Write();
00200   sigmaMassVsEta_->Write();
00201   sigmaMassVsEtaMinusErr_->Write();
00202 
00203   sigmaMassOverMassVsEtaPlusErr_->Write();
00204   sigmaMassOverMassVsEta_->Write();
00205   sigmaMassOverMassVsEtaMinusErr_->Write();
00206 
00207   canvas->Write();
00208 }
00209 
00210 // ------------ method called to for each event  ------------
00211 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00212 {
00213 }
00214 
00215 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
00216                                                  const lorentzVector& mu2,
00217                                                  const std::vector< double >& parval,
00218                                                  const double& sigmaPt1,
00219                                                  const double& sigmaPt2)
00220 {
00221   double * p = new double[(int)(parval.size())];
00222   std::vector<double>::const_iterator it = parval.begin();
00223   int id = 0;
00224   for ( ; it!=parval.end(); ++it, ++id) {
00225     p[id] = *it;
00226   }
00227   double massRes = massResolution (mu1, mu2, p, sigmaPt1, sigmaPt2);
00228   delete[] p;
00229   return massRes;
00230 }
00231 
00232 double ErrorsPropagationAnalyzer::massResolution( const lorentzVector& mu1,
00233                                                   const lorentzVector& mu2,
00234                                                   double* parval,
00235                                                   const double & sigmaPt1,
00236                                                   const double & sigmaPt2)
00237 {
00238   double mass   = (mu1+mu2).mass();
00239   double pt1    = mu1.Pt();
00240   double phi1   = mu1.Phi();
00241   double eta1   = mu1.Eta();
00242   double theta1 = 2*atan(exp(-eta1));
00243   double pt2    = mu2.Pt();
00244   double phi2   = mu2.Phi();
00245   double eta2   = mu2.Eta();
00246   double theta2 = 2*atan(exp(-eta2));
00247 
00248   // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
00249   // -----------------------------------------------------------------------------------------------------------
00250   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))-
00251                     pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
00252   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))-
00253                     pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
00254   double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
00255   double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
00256   double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
00257                        sqrt((std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)/(std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)) -
00258                        pt1*pt2*cos(theta2)/sin(theta2))/mass;
00259   double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
00260                        sqrt((std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)/(std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)) -
00261                        pt2*pt1*cos(theta1)/sin(theta1))/mass;
00262 
00263   // Resolution parameters:
00264   // ----------------------
00265   // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
00266   // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
00267   double sigma_pt1 = sigmaPt1;
00268   double sigma_pt2 = sigmaPt2;
00269   double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt1,eta1,parval );
00270   double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt2,eta2,parval );
00271   double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt1,eta1,parval );
00272   double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt2,eta2,parval );
00273   double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2( pt1, eta1, pt2, eta2, parval );
00274 
00275   // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
00276   // multiply it by pt.
00277   double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2)+
00278                          std::pow(dmdphi1*sigma_phi1,2)+std::pow(dmdphi2*sigma_phi2,2)+
00279                          std::pow(dmdcotgth1*sigma_cotgth1,2)+std::pow(dmdcotgth2*sigma_cotgth2,2)+
00280                          2*dmdpt1*dmdpt2*cov_pt1pt2);
00281 
00282   return mass_res;
00283 }
00284 
00285 
00286 void ErrorsPropagationAnalyzer::fillHistograms()
00287 {
00288   std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
00289   RootTreeHandler rootTreeHandler;
00290 
00291   typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
00292   MuonPairVector savedPair;
00293   std::vector<std::pair<int, int> > evtRun;
00294   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
00295   // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
00296 
00297   resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec = resolutionFunctionVecService( resolFitType_ );
00298   MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType_ );
00299   MuScleFitUtils::debugMassResol_ = false;
00300   MuScleFitUtils::debug = 0;
00301   MuScleFitUtils::resfind = std::vector<int>(6, 0);
00302 
00303   SigmaPt sigmaPt( parameters_, errors_ );
00304   SigmaPtDiff sigmaPtDiff;
00305 
00306   // Loop on all the pairs
00307   unsigned int i = 0;
00308   MuonPairVector::iterator it = savedPair.begin();
00309   std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
00310   for( ; it != savedPair.end(); ++it, ++i ) {
00311     double pt1 = it->first.pt();
00312     double eta1 = it->first.eta();
00313     double pt2 = it->second.pt();
00314     double eta2 = it->second.eta();
00315 
00316     if( debug_ ) {
00317       std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
00318     }
00319     // double fabsEta1 = fabs(eta1);
00320     // double fabsEta2 = fabs(eta2);
00321 
00322     if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;
00323 
00324 
00325     // double sigmaPt1 = sigmaPt( eta1 );
00326     // double sigmaPt2 = sigmaPt( eta2 );
00327 
00328     // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
00329     // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
00330     // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
00331     // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
00332 
00333     double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
00334     double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
00335     double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_);
00336     double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
00337     double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_ );
00338     double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
00339 
00340     // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
00341     // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
00342     // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
00343     double sigmaMass = massResolution( it->first, it->second, parameters_, sigmaPt1, sigmaPt2 );
00344     double sigmaMassPlusErr = massResolution( it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2 );
00345     double sigmaMassMinusErr = massResolution( it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2 );
00346 
00347     double mass = (it->first + it->second).mass();
00348     
00349     if( debug_ ) {
00350       std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
00351       std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
00352       std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
00353     }
00354 
00355     // Protections from nans
00356     if( pt1 != pt1 ) continue;
00357     if( pt2 != pt2 ) continue;
00358     if( eta1 != eta1 ) continue;
00359     if( eta2 != eta2 ) continue;
00360     if( sigmaPt1 != sigmaPt1 ) continue;
00361     if( sigmaPt2 != sigmaPt2 ) continue;
00362     if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
00363     if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
00364     if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
00365     if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
00366     if( sigmaMass != sigmaMass ) continue;
00367     if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
00368     if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
00369     if( mass == 0 ) continue;
00370 
00371     std::cout << "Muon pair number " << i << std::endl;
00372 
00373     // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
00374     // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
00375     // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
00376 
00377 
00378     if( (pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_) ) {
00379       sigmaPtVsPt_->Fill(pt1, sigmaPt1);
00380       sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
00381       sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
00382       sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
00383       sigmaMassVsPt_->Fill(pt1, sigmaMass*mass);
00384       sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr*mass);
00385       sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr*mass);
00386       sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
00387       sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
00388       sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
00389 
00390       sigmaPtVsEta_->Fill(eta1, sigmaPt1);
00391       sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
00392       sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
00393       sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
00394       sigmaMassVsEta_->Fill(eta1, sigmaMass*mass);
00395       sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr*mass);
00396       sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr*mass);
00397       sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
00398       sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
00399       sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
00400     }
00401     if( (pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_) ) {
00402       sigmaPtVsPt_->Fill(pt2, sigmaPt2);
00403       sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
00404       sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
00405       sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
00406       sigmaMassVsPt_->Fill(pt2, sigmaMass*mass);
00407       sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr*mass);
00408       sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr*mass);
00409       sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
00410       sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
00411       sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
00412 
00413       sigmaPtVsEta_->Fill(eta2, sigmaPt2);
00414       sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
00415       sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
00416       sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
00417       sigmaMassVsEta_->Fill(eta2, sigmaMass*mass);
00418       sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr*mass);
00419       sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr*mass);
00420       sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
00421       sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
00422       sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
00423     }
00424   }
00425 }
00426 
00427 //define this as a plug-in
00428 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
00429 
00430 #endif
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/MuonAnalysis/MomentumScaleCalibration/plugins/ErrorsPropagationAnalyzer.cc
