Reference Manual

Go to the documentation of this file.
00001 #ifndef RESOLUTIONANALYZER_CC
00002 #define RESOLUTIONANALYZER_CC
00003 
00004 #include "ResolutionAnalyzer.h"
00005 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
00006 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
00007 
00008 //
00009 // constants, enums and typedefs
00010 //
00011 
00012 //
00013 // static data member definitions
00014 //
00015 
00016 //
00017 // constructors and destructor
00018 //
00019 ResolutionAnalyzer::ResolutionAnalyzer(const edm::ParameterSet& iConfig) :
00020   theMuonLabel_( iConfig.getParameter<edm::InputTag>( "MuonLabel" ) ),
00021   theMuonType_( iConfig.getParameter<int>( "MuonType" ) ),
00022   theRootFileName_( iConfig.getUntrackedParameter<std::string>("OutputFileName") ),
00023   theCovariancesRootFileName_( iConfig.getUntrackedParameter<std::string>("InputFileName") ),
00024   debug_( iConfig.getUntrackedParameter<bool>( "Debug" ) ),
00025   resonance_( iConfig.getParameter<bool>("Resonance") ),
00026   readCovariances_( iConfig.getParameter<bool>( "ReadCovariances" ) ),
00027   treeFileName_( iConfig.getParameter<std::string>("InputTreeName") ),
00028   maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
00029   ptMax_( iConfig.getParameter<double>("PtMax") )
00030 {
00031   //now do what ever initialization is needed
00032 
00033   // Initial parameters values
00034   // -------------------------
00035   int resolFitType = iConfig.getParameter<int>("ResolFitType");
00036   MuScleFitUtils::ResolFitType = resolFitType;
00037   // MuScleFitUtils::resolutionFunction = resolutionFunctionArray[resolFitType];
00038   MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType );
00039   // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType];
00040   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService( resolFitType );
00041 
00042   MuScleFitUtils::parResol = iConfig.getParameter<std::vector<double> >("parResol");
00043 
00044   MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("ResFind");
00045 
00046   outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
00047   outputFile_->cd();
00048   fillHistoMap();
00049 
00050   eventCounter_ = 0;
00051 }
00052 
00053 
00054 ResolutionAnalyzer::~ResolutionAnalyzer()
00055 {
00056   outputFile_->cd();
00057   writeHistoMap();
00058   outputFile_->Close();
00059   std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
00060 }
00061 
00062 
00063 //
00064 // member functions
00065 //
00066 
00067 // ------------ method called to for each event  ------------
00068 void ResolutionAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00069   using namespace edm;
00070 
00071   std::cout << "starting" << std::endl;
00072 
00073   lorentzVector nullLorentzVector(0, 0, 0, 0);
00074 
00075   RootTreeHandler rootTreeHandler;
00076   typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
00077   MuonPairVector savedPairVector;
00078   MuonPairVector genPairVector;
00079   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &genPairVector);
00080   MuonPairVector::iterator savedPair = savedPairVector.begin();
00081   MuonPairVector::iterator genPair = genPairVector.begin();
00082   std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
00083   for( ; savedPair != savedPairVector.end(); ++savedPair, ++genPair ) {
00084 
00085   ++eventCounter_;
00086 
00087   if( (eventCounter_ % 10000) == 0 ) {
00088     std::cout << "event = " << eventCounter_ << std::endl;
00089   }
00090   
00091   lorentzVector recMu1( savedPair->first );
00092   lorentzVector recMu2( savedPair->second );
00093   
00094   if ( resonance_ ) {
00095 
00096     // Histograms with genParticles characteristics
00097     // --------------------------------------------
00098 
00099     reco::Particle::LorentzVector genMother( genPair->first + genPair->second );
00100   
00101     mapHisto_["GenMother"]->Fill( genMother );
00102     mapHisto_["DeltaGenMotherMuons"]->Fill( genPair->first, genPair->second );
00103     mapHisto_["GenMotherMuons"]->Fill( genPair->first );
00104     mapHisto_["GenMotherMuons"]->Fill( genPair->second );
00105   
00106     // Match the reco muons with the gen and sim tracks
00107     // ------------------------------------------------
00108     if(checkDeltaR(genPair->first,recMu1)){
00109       mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Pt()+recMu1.Pt())/genPair->first.Pt(),-1);
00110       mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Theta()+recMu1.Theta()),-1);
00111       mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu1,(-cos(genPair->first.Theta())/sin(genPair->first.Theta())
00112                                                             +cos(recMu1.Theta())/sin(recMu1.Theta())),-1);
00113       mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Eta()+recMu1.Eta()),-1);
00114       // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
00115       mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()),-1);
00116       recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
00117       deltaPtOverPt_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
00118       if( fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2 ) {
00119         recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
00120         deltaPtOverPtForEta12_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
00121       }
00122     }
00123     if(checkDeltaR(genPair->second,recMu2)){
00124       mapHisto_["PtResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Pt()+recMu2.Pt())/genPair->second.Pt(),+1);
00125       mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Theta()+recMu2.Theta()),+1);
00126       mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu2,(-cos(genPair->second.Theta())/sin(genPair->second.Theta())
00127                                                             +cos(recMu2.Theta())/sin(recMu2.Theta())),+1);
00128       mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Eta()+recMu2.Eta()),+1);
00129       // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
00130       mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()),+1);
00131       recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
00132       deltaPtOverPt_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
00133       if( fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2 ) {
00134         recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
00135         deltaPtOverPtForEta12_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
00136       }
00137     }  
00138 
00139     // Fill the mass resolution histograms
00140     // -----------------------------------
00141     // check if the recoMuons match the genMuons
00142     // if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
00143     if( genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
00144         checkDeltaR(genPair->first,recMu1) && checkDeltaR(genPair->second,recMu2) ) {
00145 
00146       double recoMass = (recMu1+recMu2).mass();
00147       double genMass = (genPair->first + genPair->second).mass();
00148       // first is always mu-, second is always mu+
00149       mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
00150   
00151       // Fill the reconstructed resonance
00152       reco::Particle::LorentzVector recoResonance( recMu1+recMu2 );
00153       mapHisto_["RecoResonance"]->Fill( recoResonance );
00154       mapHisto_["DeltaRecoResonanceMuons"]->Fill( recMu1, recMu2 );
00155       mapHisto_["RecoResonanceMuons"]->Fill( recMu1 );
00156       mapHisto_["RecoResonanceMuons"]->Fill( recMu2 );
00157   
00158       // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
00159       if( genMass != 0 ) {
00160         // double diffMass = (recoMass - genMass)/genMass;
00161         double diffMass = recoMass - genMass;
00162         // Fill if for both muons
00163         double pt1 = recMu1.pt();
00164         double eta1 = recMu1.eta();
00165         double pt2 = recMu2.pt();
00166         double eta2 = recMu2.eta();
00167         // This is to avoid nan
00168         if( diffMass == diffMass ) {
00169           massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
00170           massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
00171         }
00172         else {
00173           std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
00174         }
00175         // Fill with mass resolution from resolution function
00176         double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
00177         // The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
00178         // We alos take the squared value, since var = sigma^2.
00179         mapHisto_["hFunctionResolMass"]->Fill( recMu1, std::pow(massRes,2), -1 );
00180         mapHisto_["hFunctionResolMass"]->Fill( recMu2, std::pow(massRes,2), +1 );
00181       }
00182   
00183       // Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
00184       mapHisto_["hFunctionResolPt"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
00185       mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
00186       mapHisto_["hFunctionResolPhi"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
00187       mapHisto_["hFunctionResolPt"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
00188       mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
00189       mapHisto_["hFunctionResolPhi"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
00190 
00191       if( readCovariances_ ) {
00192         // Compute mass error terms
00193         // ------------------------
00194         double mass   = (recMu1+recMu2).mass();
00195         double pt1    = recMu1.Pt();
00196         double phi1   = recMu1.Phi();
00197         double eta1   = recMu1.Eta();
00198         double theta1 = 2*atan(exp(-eta1));
00199         double pt2    = recMu2.Pt();
00200         double phi2   = recMu2.Phi();
00201         double eta2   = recMu2.Eta();
00202         double theta2 = 2*atan(exp(-eta2));
00203         // Derivatives
00204         double mMu2 = MuScleFitUtils::mMu2;
00205         double dmdpt1  = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))- 
00206                           pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
00207         double dmdpt2  = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))- 
00208                           pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
00209         double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
00210         double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
00211         double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
00212                              sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) - 
00213                              pt1*pt2*cos(theta2)/sin(theta2))/mass;
00214         double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
00215                              sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) - 
00216                              pt2*pt1*cos(theta1)/sin(theta1))/mass;
00217 
00218         // Multiplied by the pt here
00219         // -------------------------
00220         double dmdpt[2] = {dmdpt1*recMu1.Pt(), dmdpt2*recMu2.Pt()};
00221         double dmdphi[2] = {dmdphi1, dmdphi2};
00222         double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
00223 
00224         // Evaluate the single terms in the mass error expression
00225 
00226         reco::Particle::LorentzVector * recMu[2] = { &recMu1, &recMu2 };
00227         int charge[2] = { -1, +1 };
00228 
00229         double fullMassRes = 0.;
00230         double massRes1 = 0.;
00231         double massRes2 = 0.;
00232         double massRes3 = 0.;
00233         double massRes4 = 0.;
00234         double massRes5 = 0.;
00235         double massRes6 = 0.;
00236         double massResPtAndPt12 = 0.;
00237 
00238         for( int i=0; i<2; ++i ) {
00239 
00240           double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
00241           double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
00242           double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
00243           double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
00244           double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
00245           double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
00246 
00247           double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
00248           double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
00249           double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
00250           double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
00251           double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
00252           double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
00253   
00254           // ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
00255           mapHisto_["MassResolutionPt"]->Fill( *(recMu[i]), ptVariance*std::pow(dmdpt[i],2), charge[i] );
00256           mapHisto_["MassResolutionCotgTheta"]->Fill( *(recMu[i]), cotgThetaVariance*std::pow(dmdcotgth[i],2), charge[i] );
00257           mapHisto_["MassResolutionPhi"]->Fill( *(recMu[i]), phiVariance*std::pow(dmdphi[i],2), charge[i] );
00258           mapHisto_["MassResolutionPt-CotgTheta"]->Fill( *(recMu[i]), 2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i], charge[i] );
00259           mapHisto_["MassResolutionPt-Phi"]->Fill( *(recMu[i]), 2*pt_phi*dmdpt[i]*dmdphi[i], charge[i] );
00260           mapHisto_["MassResolutionCotgTheta-Phi"]->Fill( *(recMu[i]), 2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i], charge[i] );
00261 
00262           mapHisto_["MassResolutionPt1-Pt2"]->Fill( *(recMu[i]), pt1_pt2*dmdpt[0]*dmdpt[1], charge[i] );
00263           mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill( *(recMu[i]), cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1], charge[i] );
00264           mapHisto_["MassResolutionPhi1-Phi2"]->Fill( *(recMu[i]), phi1_phi2*dmdphi[0]*dmdphi[1], charge[i] );
00265           // This must be filled for both configurations: 12 and 21
00266           mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1], charge[i] );
00267           mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0], charge[i] );
00268           mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[0]*dmdphi[1], charge[i] );
00269           mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[1]*dmdphi[0], charge[i] );
00270           mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1], charge[i] );
00271           mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0], charge[i] );
00272 
00273           // Sigmas for comparison
00274           mapHisto_["sigmaPtFromVariance"]->Fill( *(recMu[i]), sqrt(ptVariance), charge[i] );
00275           mapHisto_["sigmaCotgThetaFromVariance"]->Fill( *(recMu[i]), sqrt(cotgThetaVariance), charge[i] );
00276           mapHisto_["sigmaPhiFromVariance"]->Fill( *(recMu[i]), sqrt(phiVariance), charge[i] );
00277 
00278           // Pt term from function
00279           mapHisto_["MassResolutionPtFromFunction"]->Fill( *(recMu[i]), ( MuScleFitUtils::resolutionFunctionForVec->sigmaPt((recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol) )*std::pow(dmdpt[i],2), charge[i] );
00280 
00281           fullMassRes +=
00282             ptVariance*std::pow(dmdpt[i],2) +
00283             cotgThetaVariance*std::pow(dmdcotgth[i],2) +
00284             phiVariance*std::pow(dmdphi[i],2) +
00285 
00286             // These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
00287             2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
00288             2*pt_phi*dmdpt[i]*dmdphi[i] +
00289             2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
00290 
00291             pt1_pt2*dmdpt[0]*dmdpt[1] +
00292             cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
00293             phi1_phi2*dmdphi[0]*dmdphi[1] +
00294 
00295             // These are filled twice, because of the two combinations
00296             pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
00297             pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
00298             pt12_phi21*dmdpt[0]*dmdphi[1] +
00299             pt12_phi21*dmdpt[1]*dmdphi[0] +
00300             cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
00301             cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
00302 
00303           massRes1 += ptVariance*std::pow(dmdpt[i],2);
00304           massRes2 += ptVariance*std::pow(dmdpt[i],2) +
00305             cotgThetaVariance*std::pow(dmdcotgth[i],2);
00306           massRes3 += ptVariance*std::pow(dmdpt[i],2) +
00307             cotgThetaVariance*std::pow(dmdcotgth[i],2) +
00308             phiVariance*std::pow(dmdphi[i],2);
00309           massRes4 += ptVariance*std::pow(dmdpt[i],2) +
00310             cotgThetaVariance*std::pow(dmdcotgth[i],2) +
00311             phiVariance*std::pow(dmdphi[i],2) +
00312             pt1_pt2*dmdpt[0]*dmdpt[1] +
00313             2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
00314             2*pt_phi*dmdpt[i]*dmdphi[i] +
00315             2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i];
00316           massRes5 += ptVariance*std::pow(dmdpt[i],2) +
00317             cotgThetaVariance*std::pow(dmdcotgth[i],2) +
00318             phiVariance*std::pow(dmdphi[i],2) +
00319             pt1_pt2*dmdpt[0]*dmdpt[1] +
00320             2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
00321             2*pt_phi*dmdpt[i]*dmdphi[i] +
00322             2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
00323             cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
00324             phi1_phi2*dmdphi[0]*dmdphi[1];
00325           massRes6 += ptVariance*std::pow(dmdpt[i],2) +
00326             cotgThetaVariance*std::pow(dmdcotgth[i],2) +
00327             phiVariance*std::pow(dmdphi[i],2) +
00328             pt1_pt2*dmdpt[0]*dmdpt[1] +
00329             2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
00330             2*pt_phi*dmdpt[i]*dmdphi[i] +
00331             2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
00332             cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
00333             phi1_phi2*dmdphi[0]*dmdphi[1] +
00334             pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
00335             pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
00336             pt12_phi21*dmdpt[0]*dmdphi[1] +
00337             pt12_phi21*dmdpt[1]*dmdphi[0] +
00338             cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
00339             cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
00340 
00341           massResPtAndPt12 += ptVariance*std::pow(dmdpt[i],2) + pt1_pt2*dmdpt[0]*dmdpt[1];
00342 
00343           // Derivatives
00344           mapHisto_["DerivativePt"]->Fill( *(recMu[i]), dmdpt[i], charge[i] );
00345           mapHisto_["DerivativeCotgTheta"]->Fill( *(recMu[i]), dmdcotgth[i], charge[i] );
00346           mapHisto_["DerivativePhi"]->Fill( *(recMu[i]), dmdphi[i], charge[i] );
00347         }
00348         // Fill the complete resolution function with covariance terms
00349         mapHisto_["FullMassResolution"]->Fill( *(recMu[0]), fullMassRes, charge[0]);
00350         mapHisto_["FullMassResolution"]->Fill( *(recMu[1]), fullMassRes, charge[1]);
00351 
00352         mapHisto_["MassRes1"]->Fill( *(recMu[0]), massRes1, charge[0] );
00353         mapHisto_["MassRes1"]->Fill( *(recMu[1]), massRes1, charge[1] );
00354         mapHisto_["MassRes2"]->Fill( *(recMu[0]), massRes2, charge[0] );
00355         mapHisto_["MassRes2"]->Fill( *(recMu[1]), massRes2, charge[1] );
00356         mapHisto_["MassRes3"]->Fill( *(recMu[0]), massRes3, charge[0] );
00357         mapHisto_["MassRes3"]->Fill( *(recMu[1]), massRes3, charge[1] );
00358         mapHisto_["MassRes4"]->Fill( *(recMu[0]), massRes4, charge[0] );
00359         mapHisto_["MassRes4"]->Fill( *(recMu[1]), massRes4, charge[1] );
00360         mapHisto_["MassRes5"]->Fill( *(recMu[0]), massRes5, charge[0] );
00361         mapHisto_["MassRes5"]->Fill( *(recMu[1]), massRes5, charge[1] );
00362         mapHisto_["MassRes6"]->Fill( *(recMu[0]), massRes6, charge[0] );
00363         mapHisto_["MassRes6"]->Fill( *(recMu[1]), massRes6, charge[1] );
00364         mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[0]), massResPtAndPt12, charge[0] );
00365         mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[1]), massResPtAndPt12, charge[1] );
00366       }
00367       else {
00368         // Fill the covariances histograms
00369         mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
00370       }
00371     }
00372   } // end if resonance
00373   }
00374 //   else {
00375 // 
00376 //     // Loop on the recMuons
00377 //     std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
00378 //     for ( ; recMuon!=muons.end(); ++recMuon ) {  
00379 //       int charge = recMuon->charge();
00380 // 
00381 //       lorentzVector recMu(recMuon->p4());
00382 // 
00383 //       // Find the matching MC muon
00384 //       const HepMC::GenEvent* Evt = evtMC->GetEvent();
00385 //       //Loop on generated particles
00386 //       std::map<double, lorentzVector> genAssocMap;
00387 //       HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
00388 //       for( ; part!=Evt->particles_end(); ++part ) {
00389 //         if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
00390 //           lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00391 //                                                (*part)->momentum().pz(),(*part)->momentum().e()));
00392 // 
00393 //           double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
00394 //                                ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
00395 // 
00396 //           // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
00397 //           // Only in this case we consider it matching.
00398 //           if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
00399 //         }
00400 //       }
00401 //       // Take the closest in deltaR
00402 //       lorentzVector genMu(genAssocMap.begin()->second);
00403 // 
00404 //       // Histograms with genParticles characteristics
00405 //       // --------------------------------------------
00406 // 
00407 //       if(checkDeltaR(genMu,recMu)){
00408 //         mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
00409 //         mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
00410 //         mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
00411 //                                                              +cos(recMu.Theta())/sin(recMu.Theta())),charge);
00412 //         mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
00413 //         mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
00414 //       }
00415 //     }
00416 }
00417 
00418 void ResolutionAnalyzer::fillHistoMap() {
00419 
00420   outputFile_->cd();
00421 
00422   // Resonances
00423   // If no Z is required, use a smaller mass range.
00424   double minMass = 0.;
00425   double maxMass = 200.;
00426   if( MuScleFitUtils::resfind[0] != 1 ) {
00427     maxMass = 30.;
00428   }
00429   mapHisto_["GenMother"]               = new HParticle(outputFile_, "GenMother", minMass, maxMass);
00430   mapHisto_["SimResonance"]            = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
00431   mapHisto_["RecoResonance"]           = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
00432   
00433   // Resonance muons
00434   mapHisto_["GenMotherMuons"]          = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
00435   mapHisto_["SimResonanceMuons"]       = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
00436   mapHisto_["RecoResonanceMuons"]      = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
00437   
00438   // Deltas between resonance muons
00439   mapHisto_["DeltaGenMotherMuons"]     = new HDelta (outputFile_, "DeltaGenMotherMuons");
00440   mapHisto_["DeltaSimResonanceMuons"]  = new HDelta (outputFile_, "DeltaSimResonanceMuons");
00441   mapHisto_["DeltaRecoResonanceMuons"] = new HDelta (outputFile_, "DeltaRecoResonanceMuons");
00442   
00443   //   //Reconstructed muon kinematics
00444   //   //-----------------------------
00445   //   mapHisto_["hRecBestMu"]             = new HParticle         ("hRecBestMu");
00446   //   mapHisto_["hRecBestMu_Acc"]         = new HParticle         ("hRecBestMu_Acc"); 
00447   //   mapHisto_["hDeltaRecBestMu"]        = new HDelta            ("hDeltaRecBestMu");
00448   
00449   //   mapHisto_["hRecBestRes"]            = new HParticle         ("hRecBestRes");
00450   //   mapHisto_["hRecBestRes_Acc"]        = new HParticle         ("hRecBestRes_Acc"); 
00451   //   mapHisto_["hRecBestResVSMu"]        = new HMassVSPart       ("hRecBestResVSMu");
00452   
00453   //Resolution VS muon kinematic
00454   //----------------------------
00455   mapHisto_["PtResolutionGenVSMu"]        = new HResolutionVSPart (outputFile_, "PtResolutionGenVSMu");
00456   mapHisto_["PtResolutionSimVSMu"]        = new HResolutionVSPart (outputFile_, "PtResolutionSimVSMu");
00457   mapHisto_["EtaResolutionGenVSMu"]       = new HResolutionVSPart (outputFile_, "EtaResolutionGenVSMu");
00458   mapHisto_["EtaResolutionSimVSMu"]       = new HResolutionVSPart (outputFile_, "EtaResolutionSimVSMu");
00459   mapHisto_["ThetaResolutionGenVSMu"]     = new HResolutionVSPart (outputFile_, "ThetaResolutionGenVSMu");
00460   mapHisto_["ThetaResolutionSimVSMu"]     = new HResolutionVSPart (outputFile_, "ThetaResolutionSimVSMu");
00461   mapHisto_["CotgThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
00462   mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionSimVSMu");
00463   mapHisto_["PhiResolutionGenVSMu"]       = new HResolutionVSPart (outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
00464   mapHisto_["PhiResolutionSimVSMu"]       = new HResolutionVSPart (outputFile_, "PhiResolutionSimVSMu");
00465 
00466   // Covariances between muons kinematic quantities
00467   // ----------------------------------------------
00468   double ptMax = ptMax_;
00469 
00470   // Mass resolution
00471   // ---------------
00472   mapHisto_["MassResolution"] = new HMassResolutionVSPart (outputFile_,"MassResolution");
00473   
00474   //  mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
00475   
00476   // Mass resolution vs (pt, eta) of the muons from MC
00477   massResolutionVsPtEta_ = new HCovarianceVSxy ( "Mass", "Mass", 100, 0., ptMax, 60, -3, 3 );
00478   // Mass resolution vs (pt, eta) of the muons from function
00479   recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
00480   recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
00481   deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
00482   deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
00483 
00484   // Muons resolutions from resolution functions
00485   // -------------------------------------------
00486   int totBinsY = 60;
00487   mapHisto_["hFunctionResolMass"]      = new HFunctionResolution (outputFile_, "hFunctionResolMass", ptMax, totBinsY);
00488   mapHisto_["hFunctionResolPt"]        = new HFunctionResolution (outputFile_, "hFunctionResolPt", ptMax, totBinsY);
00489   mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution (outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
00490   mapHisto_["hFunctionResolPhi"]       = new HFunctionResolution (outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
00491 
00492   if( readCovariances_ ) {
00493     // Covariances read from file. Used to compare the terms in the expression of mass error
00494     mapHisto_["ReadCovariances"] = new HCovarianceVSParts ( theCovariancesRootFileName_, "Covariance" );
00495 
00496     // Variances
00497     mapHisto_["MassResolutionPt"]                    = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPt", ptMax);
00498     mapHisto_["MassResolutionCotgTheta"]             = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassCotgTheta", ptMax);
00499     mapHisto_["MassResolutionPhi"]                   = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPhi", ptMax);
00500     // Covariances
00501     mapHisto_["MassResolutionPt-CotgTheta"]          = new HFunctionResolution(outputFile_,"functionResolMassPt-CotgTheta", ptMax, totBinsY);
00502     mapHisto_["MassResolutionPt-Phi"]                = new HFunctionResolution(outputFile_,"functionResolMassPt-Phi", ptMax, totBinsY);
00503     mapHisto_["MassResolutionCotgTheta-Phi"]         = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta-Phi", ptMax, totBinsY);
00504     mapHisto_["MassResolutionPt1-Pt2"]               = new HFunctionResolution(outputFile_,"functionResolMassPt1-Pt2", ptMax, totBinsY);
00505     mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
00506     mapHisto_["MassResolutionPhi1-Phi2"]             = new HFunctionResolution(outputFile_,"functionResolMassPhi1-Phi2", ptMax, totBinsY);
00507     mapHisto_["MassResolutionPt12-CotgTheta21"]      = new HFunctionResolution(outputFile_,"functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
00508     mapHisto_["MassResolutionPt12-Phi21"]            = new HFunctionResolution(outputFile_,"functionResolMassPt12-Phi21", ptMax, totBinsY);
00509     mapHisto_["MassResolutionCotgTheta12-Phi21"]     = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
00510 
00511     mapHisto_["sigmaPtFromVariance"]                 = new HFunctionResolution(outputFile_,"sigmaPtFromVariance", ptMax, totBinsY);
00512     mapHisto_["sigmaCotgThetaFromVariance"]          = new HFunctionResolution(outputFile_,"sigmaCotgThetaFromVariance", ptMax, totBinsY);
00513     mapHisto_["sigmaPhiFromVariance"]                = new HFunctionResolution(outputFile_,"sigmaPhiFromVariance", ptMax, totBinsY);
00514 
00515     // Derivatives
00516     mapHisto_["DerivativePt"]                        = new HFunctionResolution(outputFile_,"derivativePt", ptMax);
00517     mapHisto_["DerivativeCotgTheta"]                 = new HFunctionResolution(outputFile_,"derivativeCotgTheta", ptMax);
00518     mapHisto_["DerivativePhi"]                       = new HFunctionResolution(outputFile_,"derivativePhi", ptMax);
00519 
00520     // Pt term from function
00521     mapHisto_["MassResolutionPtFromFunction"]        = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPtFromFunction", ptMax);
00522 
00523     mapHisto_["FullMassResolution"]                  = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
00524     mapHisto_["MassRes1"]                            = new HFunctionResolution(outputFile_, "massRes1", ptMax);
00525     mapHisto_["MassRes2"]                            = new HFunctionResolution(outputFile_, "massRes2", ptMax);
00526     mapHisto_["MassRes3"]                            = new HFunctionResolution(outputFile_, "massRes3", ptMax);
00527     mapHisto_["MassRes4"]                            = new HFunctionResolution(outputFile_, "massRes4", ptMax);
00528     mapHisto_["MassRes5"]                            = new HFunctionResolution(outputFile_, "massRes5", ptMax);
00529     mapHisto_["MassRes6"]                            = new HFunctionResolution(outputFile_, "massRes6", ptMax);
00530     mapHisto_["MassResPtAndPt12"]                    = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
00531   }
00532   else {
00533     mapHisto_["Covariances"] = new HCovarianceVSParts ( outputFile_, "Covariance", ptMax );
00534   }
00535 }
00536 
00537 void ResolutionAnalyzer::writeHistoMap() {
00538   for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin(); 
00539        histo!=mapHisto_.end(); histo++) {
00540     (*histo).second->Write();
00541   }
00542   outputFile_->cd();
00543   massResolutionVsPtEta_->Write();
00544   recoPtVsgenPt_->Write();
00545   recoPtVsgenPtEta12_->Write();
00546   deltaPtOverPt_->Write();
00547   deltaPtOverPtForEta12_->Write();
00548 }
00549 
00550 bool ResolutionAnalyzer::checkDeltaR(const reco::Particle::LorentzVector & genMu, const reco::Particle::LorentzVector & recMu){
00551   //first is always mu-, second is always mu+
00552   double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) +
00553                        ((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
00554   if(deltaR<0.01)
00555     return true;
00556   else if (debug_>0)
00557     std::cout<<"Reco muon "<<recMu<<" with eta "<<recMu.Eta()<<" and phi "<<recMu.Phi()<<std::endl
00558              <<" DOES NOT MATCH with generated muon from resonance: "<<std::endl
00559              <<genMu<<" with eta "<<genMu.Eta()<<" and phi "<<genMu.Phi()<<std::endl;
00560   return false;
00561 }
00562 
00563 //define this as a plug-in
00564 DEFINE_FWK_MODULE(ResolutionAnalyzer);
00565 
00566 #endif // RESOLUTIONANALYZER_CC
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/MuonAnalysis/MomentumScaleCalibration/plugins/ResolutionAnalyzer.cc
