|
|
|
#include <MuonAnalysis/MomentumScaleCalibration/plugins/ResolutionAnalyzer.cc>
Public Member Functions | |
| ResolutionAnalyzer (const edm::ParameterSet &) | |
| ~ResolutionAnalyzer () | |
Private Member Functions | |
| virtual void | analyze (const edm::Event &, const edm::EventSetup &) |
| bool | checkDeltaR (const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recMu) |
| Returns true if the two particles have DeltaR < cut. | |
| virtual void | endJob () |
| void | fillHistoMap () |
| Used to fill the map with the histograms needed. | |
| template<typename T > | |
| std::vector< reco::LeafCandidate > | fillMuonCollection (const std::vector< T > &tracks) |
| void | writeHistoMap () |
| Writes the histograms in the map. | |
Private Attributes | |
| bool | debug_ |
| TH1D * | deltaPtOverPt_ |
| TH1D * | deltaPtOverPtForEta12_ |
| int | eventCounter_ |
| std::map< std::string, Histograms * > | mapHisto_ |
| HCovarianceVSxy * | massResolutionVsPtEta_ |
| int32_t | maxEvents_ |
| TFile * | outputFile_ |
| double | ptMax_ |
| bool | readCovariances_ |
| TH2D * | recoPtVsgenPt_ |
| TH2D * | recoPtVsgenPtEta12_ |
| bool | resonance_ |
| std::string | theCovariancesRootFileName_ |
| edm::InputTag | theMuonLabel_ |
| int | theMuonType_ |
| std::string | theRootFileName_ |
| TString | treeFileName_ |
Description: <one line="" class="" summary>="">
Implementation: <Notes on="" implementation>="">
Definition at line 66 of file ResolutionAnalyzer.h.
| ResolutionAnalyzer::ResolutionAnalyzer | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 19 of file ResolutionAnalyzer.cc.
References eventCounter_, fillHistoMap(), edm::ParameterSet::getParameter(), outputFile_, MuScleFitUtils::parResol, MuScleFitUtils::resfind, MuScleFitUtils::ResolFitType, MuScleFitUtils::resolutionFunction, MuScleFitUtils::resolutionFunctionForVec, resolutionFunctionService(), resolutionFunctionVecService(), and theRootFileName_.
: theMuonLabel_( iConfig.getParameter<edm::InputTag>( "MuonLabel" ) ), theMuonType_( iConfig.getParameter<int>( "MuonType" ) ), theRootFileName_( iConfig.getUntrackedParameter<std::string>("OutputFileName") ), theCovariancesRootFileName_( iConfig.getUntrackedParameter<std::string>("InputFileName") ), debug_( iConfig.getUntrackedParameter<bool>( "Debug" ) ), resonance_( iConfig.getParameter<bool>("Resonance") ), readCovariances_( iConfig.getParameter<bool>( "ReadCovariances" ) ), treeFileName_( iConfig.getParameter<std::string>("InputTreeName") ), maxEvents_( iConfig.getParameter<int>("MaxEvents") ), ptMax_( iConfig.getParameter<double>("PtMax") ) { //now do what ever initialization is needed // Initial parameters values // ------------------------- int resolFitType = iConfig.getParameter<int>("ResolFitType"); MuScleFitUtils::ResolFitType = resolFitType; // MuScleFitUtils::resolutionFunction = resolutionFunctionArray[resolFitType]; MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType ); // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType]; MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService( resolFitType ); MuScleFitUtils::parResol = iConfig.getParameter<std::vector<double> >("parResol"); MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("ResFind"); outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE"); outputFile_->cd(); fillHistoMap(); eventCounter_ = 0; }
| ResolutionAnalyzer::~ResolutionAnalyzer | ( | ) |
Definition at line 54 of file ResolutionAnalyzer.cc.
References gather_cfg::cout, eventCounter_, outputFile_, and writeHistoMap().
{
outputFile_->cd();
writeHistoMap();
outputFile_->Close();
std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
}
| void ResolutionAnalyzer::analyze | ( | const edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [private, virtual] |
Implements edm::EDAnalyzer.
Definition at line 68 of file ResolutionAnalyzer.cc.
References DeDxDiscriminatorTools::charge(), checkDeltaR(), funct::cos(), gather_cfg::cout, MuScleFitUtils::deltaPhiNoFabs(), deltaPtOverPt_, deltaPtOverPtForEta12_, reco::tau::disc::Eta(), eventCounter_, create_public_lumi_plots::exp, HCovarianceVSxy::Fill(), i, mapHisto_, scaleCards::mass, MuScleFitUtils::massResolution(), massResolutionVsPtEta_, maxEvents_, MuScleFitUtils::mMu2, MuScleFitUtils::parResol, funct::pow(), reco::tau::disc::Pt(), readCovariances_, recoPtVsgenPt_, recoPtVsgenPtEta12_, MuScleFitUtils::resolutionFunctionForVec, resonance_, funct::sin(), mathSSE::sqrt(), and treeFileName_.
{
using namespace edm;
std::cout << "starting" << std::endl;
lorentzVector nullLorentzVector(0, 0, 0, 0);
RootTreeHandler rootTreeHandler;
typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
MuonPairVector savedPairVector;
MuonPairVector genPairVector;
std::vector<std::pair<int, int> > evtRun;
rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &evtRun, &genPairVector);
MuonPairVector::iterator savedPair = savedPairVector.begin();
MuonPairVector::iterator genPair = genPairVector.begin();
std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
for( ; savedPair != savedPairVector.end(); ++savedPair, ++genPair ) {
++eventCounter_;
if( (eventCounter_ % 10000) == 0 ) {
std::cout << "event = " << eventCounter_ << std::endl;
}
lorentzVector recMu1( savedPair->first );
lorentzVector recMu2( savedPair->second );
if ( resonance_ ) {
// Histograms with genParticles characteristics
// --------------------------------------------
reco::Particle::LorentzVector genMother( genPair->first + genPair->second );
mapHisto_["GenMother"]->Fill( genMother );
mapHisto_["DeltaGenMotherMuons"]->Fill( genPair->first, genPair->second );
mapHisto_["GenMotherMuons"]->Fill( genPair->first );
mapHisto_["GenMotherMuons"]->Fill( genPair->second );
// Match the reco muons with the gen and sim tracks
// ------------------------------------------------
if(checkDeltaR(genPair->first,recMu1)){
mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Pt()+recMu1.Pt())/genPair->first.Pt(),-1);
mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Theta()+recMu1.Theta()),-1);
mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu1,(-cos(genPair->first.Theta())/sin(genPair->first.Theta())
+cos(recMu1.Theta())/sin(recMu1.Theta())),-1);
mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Eta()+recMu1.Eta()),-1);
// mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()),-1);
recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
deltaPtOverPt_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
if( fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2 ) {
recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
deltaPtOverPtForEta12_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
}
}
if(checkDeltaR(genPair->second,recMu2)){
mapHisto_["PtResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Pt()+recMu2.Pt())/genPair->second.Pt(),+1);
mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Theta()+recMu2.Theta()),+1);
mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu2,(-cos(genPair->second.Theta())/sin(genPair->second.Theta())
+cos(recMu2.Theta())/sin(recMu2.Theta())),+1);
mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Eta()+recMu2.Eta()),+1);
// mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()),+1);
recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
deltaPtOverPt_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
if( fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2 ) {
recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
deltaPtOverPtForEta12_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
}
}
// Fill the mass resolution histograms
// -----------------------------------
// check if the recoMuons match the genMuons
// if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
if( genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
checkDeltaR(genPair->first,recMu1) && checkDeltaR(genPair->second,recMu2) ) {
double recoMass = (recMu1+recMu2).mass();
double genMass = (genPair->first + genPair->second).mass();
// first is always mu-, second is always mu+
mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
// Fill the reconstructed resonance
reco::Particle::LorentzVector recoResonance( recMu1+recMu2 );
mapHisto_["RecoResonance"]->Fill( recoResonance );
mapHisto_["DeltaRecoResonanceMuons"]->Fill( recMu1, recMu2 );
mapHisto_["RecoResonanceMuons"]->Fill( recMu1 );
mapHisto_["RecoResonanceMuons"]->Fill( recMu2 );
// Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
if( genMass != 0 ) {
// double diffMass = (recoMass - genMass)/genMass;
double diffMass = recoMass - genMass;
// Fill if for both muons
double pt1 = recMu1.pt();
double eta1 = recMu1.eta();
double pt2 = recMu2.pt();
double eta2 = recMu2.eta();
// This is to avoid nan
if( diffMass == diffMass ) {
massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
}
else {
std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
}
// Fill with mass resolution from resolution function
double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
// The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
// We alos take the squared value, since var = sigma^2.
mapHisto_["hFunctionResolMass"]->Fill( recMu1, std::pow(massRes,2), -1 );
mapHisto_["hFunctionResolMass"]->Fill( recMu2, std::pow(massRes,2), +1 );
}
// Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
mapHisto_["hFunctionResolPt"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
mapHisto_["hFunctionResolPhi"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
mapHisto_["hFunctionResolPt"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
mapHisto_["hFunctionResolPhi"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
if( readCovariances_ ) {
// Compute mass error terms
// ------------------------
double mass = (recMu1+recMu2).mass();
double pt1 = recMu1.Pt();
double phi1 = recMu1.Phi();
double eta1 = recMu1.Eta();
double theta1 = 2*atan(exp(-eta1));
double pt2 = recMu2.Pt();
double phi2 = recMu2.Phi();
double eta2 = recMu2.Eta();
double theta2 = 2*atan(exp(-eta2));
// Derivatives
double mMu2 = MuScleFitUtils::mMu2;
double dmdpt1 = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))-
pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
double dmdpt2 = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))-
pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) -
pt1*pt2*cos(theta2)/sin(theta2))/mass;
double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) -
pt2*pt1*cos(theta1)/sin(theta1))/mass;
// Multiplied by the pt here
// -------------------------
double dmdpt[2] = {dmdpt1*recMu1.Pt(), dmdpt2*recMu2.Pt()};
double dmdphi[2] = {dmdphi1, dmdphi2};
double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
// Evaluate the single terms in the mass error expression
reco::Particle::LorentzVector * recMu[2] = { &recMu1, &recMu2 };
int charge[2] = { -1, +1 };
double fullMassRes = 0.;
double massRes1 = 0.;
double massRes2 = 0.;
double massRes3 = 0.;
double massRes4 = 0.;
double massRes5 = 0.;
double massRes6 = 0.;
double massResPtAndPt12 = 0.;
for( int i=0; i<2; ++i ) {
double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
// ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
mapHisto_["MassResolutionPt"]->Fill( *(recMu[i]), ptVariance*std::pow(dmdpt[i],2), charge[i] );
mapHisto_["MassResolutionCotgTheta"]->Fill( *(recMu[i]), cotgThetaVariance*std::pow(dmdcotgth[i],2), charge[i] );
mapHisto_["MassResolutionPhi"]->Fill( *(recMu[i]), phiVariance*std::pow(dmdphi[i],2), charge[i] );
mapHisto_["MassResolutionPt-CotgTheta"]->Fill( *(recMu[i]), 2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i], charge[i] );
mapHisto_["MassResolutionPt-Phi"]->Fill( *(recMu[i]), 2*pt_phi*dmdpt[i]*dmdphi[i], charge[i] );
mapHisto_["MassResolutionCotgTheta-Phi"]->Fill( *(recMu[i]), 2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i], charge[i] );
mapHisto_["MassResolutionPt1-Pt2"]->Fill( *(recMu[i]), pt1_pt2*dmdpt[0]*dmdpt[1], charge[i] );
mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill( *(recMu[i]), cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1], charge[i] );
mapHisto_["MassResolutionPhi1-Phi2"]->Fill( *(recMu[i]), phi1_phi2*dmdphi[0]*dmdphi[1], charge[i] );
// This must be filled for both configurations: 12 and 21
mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1], charge[i] );
mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0], charge[i] );
mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[0]*dmdphi[1], charge[i] );
mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[1]*dmdphi[0], charge[i] );
mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1], charge[i] );
mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0], charge[i] );
// Sigmas for comparison
mapHisto_["sigmaPtFromVariance"]->Fill( *(recMu[i]), sqrt(ptVariance), charge[i] );
mapHisto_["sigmaCotgThetaFromVariance"]->Fill( *(recMu[i]), sqrt(cotgThetaVariance), charge[i] );
mapHisto_["sigmaPhiFromVariance"]->Fill( *(recMu[i]), sqrt(phiVariance), charge[i] );
// Pt term from function
mapHisto_["MassResolutionPtFromFunction"]->Fill( *(recMu[i]), ( MuScleFitUtils::resolutionFunctionForVec->sigmaPt((recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol) )*std::pow(dmdpt[i],2), charge[i] );
fullMassRes +=
ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2) +
phiVariance*std::pow(dmdphi[i],2) +
// These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
2*pt_phi*dmdpt[i]*dmdphi[i] +
2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
pt1_pt2*dmdpt[0]*dmdpt[1] +
cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
phi1_phi2*dmdphi[0]*dmdphi[1] +
// These are filled twice, because of the two combinations
pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
pt12_phi21*dmdpt[0]*dmdphi[1] +
pt12_phi21*dmdpt[1]*dmdphi[0] +
cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
massRes1 += ptVariance*std::pow(dmdpt[i],2);
massRes2 += ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2);
massRes3 += ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2) +
phiVariance*std::pow(dmdphi[i],2);
massRes4 += ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2) +
phiVariance*std::pow(dmdphi[i],2) +
pt1_pt2*dmdpt[0]*dmdpt[1] +
2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
2*pt_phi*dmdpt[i]*dmdphi[i] +
2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i];
massRes5 += ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2) +
phiVariance*std::pow(dmdphi[i],2) +
pt1_pt2*dmdpt[0]*dmdpt[1] +
2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
2*pt_phi*dmdpt[i]*dmdphi[i] +
2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
phi1_phi2*dmdphi[0]*dmdphi[1];
massRes6 += ptVariance*std::pow(dmdpt[i],2) +
cotgThetaVariance*std::pow(dmdcotgth[i],2) +
phiVariance*std::pow(dmdphi[i],2) +
pt1_pt2*dmdpt[0]*dmdpt[1] +
2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
2*pt_phi*dmdpt[i]*dmdphi[i] +
2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
phi1_phi2*dmdphi[0]*dmdphi[1] +
pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
pt12_phi21*dmdpt[0]*dmdphi[1] +
pt12_phi21*dmdpt[1]*dmdphi[0] +
cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
massResPtAndPt12 += ptVariance*std::pow(dmdpt[i],2) + pt1_pt2*dmdpt[0]*dmdpt[1];
// Derivatives
mapHisto_["DerivativePt"]->Fill( *(recMu[i]), dmdpt[i], charge[i] );
mapHisto_["DerivativeCotgTheta"]->Fill( *(recMu[i]), dmdcotgth[i], charge[i] );
mapHisto_["DerivativePhi"]->Fill( *(recMu[i]), dmdphi[i], charge[i] );
}
// Fill the complete resolution function with covariance terms
mapHisto_["FullMassResolution"]->Fill( *(recMu[0]), fullMassRes, charge[0]);
mapHisto_["FullMassResolution"]->Fill( *(recMu[1]), fullMassRes, charge[1]);
mapHisto_["MassRes1"]->Fill( *(recMu[0]), massRes1, charge[0] );
mapHisto_["MassRes1"]->Fill( *(recMu[1]), massRes1, charge[1] );
mapHisto_["MassRes2"]->Fill( *(recMu[0]), massRes2, charge[0] );
mapHisto_["MassRes2"]->Fill( *(recMu[1]), massRes2, charge[1] );
mapHisto_["MassRes3"]->Fill( *(recMu[0]), massRes3, charge[0] );
mapHisto_["MassRes3"]->Fill( *(recMu[1]), massRes3, charge[1] );
mapHisto_["MassRes4"]->Fill( *(recMu[0]), massRes4, charge[0] );
mapHisto_["MassRes4"]->Fill( *(recMu[1]), massRes4, charge[1] );
mapHisto_["MassRes5"]->Fill( *(recMu[0]), massRes5, charge[0] );
mapHisto_["MassRes5"]->Fill( *(recMu[1]), massRes5, charge[1] );
mapHisto_["MassRes6"]->Fill( *(recMu[0]), massRes6, charge[0] );
mapHisto_["MassRes6"]->Fill( *(recMu[1]), massRes6, charge[1] );
mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[0]), massResPtAndPt12, charge[0] );
mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[1]), massResPtAndPt12, charge[1] );
}
else {
// Fill the covariances histograms
mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
}
}
} // end if resonance
}
// else {
//
// // Loop on the recMuons
// std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
// for ( ; recMuon!=muons.end(); ++recMuon ) {
// int charge = recMuon->charge();
//
// lorentzVector recMu(recMuon->p4());
//
// // Find the matching MC muon
// const HepMC::GenEvent* Evt = evtMC->GetEvent();
// //Loop on generated particles
// std::map<double, lorentzVector> genAssocMap;
// HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
// for( ; part!=Evt->particles_end(); ++part ) {
// if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
// lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
// (*part)->momentum().pz(),(*part)->momentum().e()));
//
// double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
// ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
//
// // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
// // Only in this case we consider it matching.
// if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
// }
// }
// // Take the closest in deltaR
// lorentzVector genMu(genAssocMap.begin()->second);
//
// // Histograms with genParticles characteristics
// // --------------------------------------------
//
// if(checkDeltaR(genMu,recMu)){
// mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
// mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
// mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
// +cos(recMu.Theta())/sin(recMu.Theta())),charge);
// mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
// mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
// }
// }
}
| bool ResolutionAnalyzer::checkDeltaR | ( | const reco::Particle::LorentzVector & | genMu, |
| const reco::Particle::LorentzVector & | recMu | ||
| ) | [private] |
Returns true if the two particles have DeltaR < cut.
Definition at line 552 of file ResolutionAnalyzer.cc.
References gather_cfg::cout, debug_, MuScleFitUtils::deltaPhi(), deltaR(), and mathSSE::sqrt().
Referenced by analyze().
{
//first is always mu-, second is always mu+
double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) +
((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
if(deltaR<0.01)
return true;
else if (debug_>0)
std::cout<<"Reco muon "<<recMu<<" with eta "<<recMu.Eta()<<" and phi "<<recMu.Phi()<<std::endl
<<" DOES NOT MATCH with generated muon from resonance: "<<std::endl
<<genMu<<" with eta "<<genMu.Eta()<<" and phi "<<genMu.Phi()<<std::endl;
return false;
}
| virtual void ResolutionAnalyzer::endJob | ( | void | ) | [inline, private, virtual] |
| void ResolutionAnalyzer::fillHistoMap | ( | ) | [private] |
Used to fill the map with the histograms needed.
Definition at line 420 of file ResolutionAnalyzer.cc.
References deltaPtOverPt_, deltaPtOverPtForEta12_, mapHisto_, massResolutionVsPtEta_, outputFile_, jptDQMConfig_cff::ptMax, ptMax_, readCovariances_, recoPtVsgenPt_, recoPtVsgenPtEta12_, MuScleFitUtils::resfind, and theCovariancesRootFileName_.
Referenced by ResolutionAnalyzer().
{
outputFile_->cd();
// Resonances
// If no Z is required, use a smaller mass range.
double minMass = 0.;
double maxMass = 200.;
if( MuScleFitUtils::resfind[0] != 1 ) {
maxMass = 30.;
}
mapHisto_["GenMother"] = new HParticle(outputFile_, "GenMother", minMass, maxMass);
mapHisto_["SimResonance"] = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
mapHisto_["RecoResonance"] = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
// Resonance muons
mapHisto_["GenMotherMuons"] = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
mapHisto_["SimResonanceMuons"] = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
mapHisto_["RecoResonanceMuons"] = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
// Deltas between resonance muons
mapHisto_["DeltaGenMotherMuons"] = new HDelta (outputFile_, "DeltaGenMotherMuons");
mapHisto_["DeltaSimResonanceMuons"] = new HDelta (outputFile_, "DeltaSimResonanceMuons");
mapHisto_["DeltaRecoResonanceMuons"] = new HDelta (outputFile_, "DeltaRecoResonanceMuons");
// //Reconstructed muon kinematics
// //-----------------------------
// mapHisto_["hRecBestMu"] = new HParticle ("hRecBestMu");
// mapHisto_["hRecBestMu_Acc"] = new HParticle ("hRecBestMu_Acc");
// mapHisto_["hDeltaRecBestMu"] = new HDelta ("hDeltaRecBestMu");
// mapHisto_["hRecBestRes"] = new HParticle ("hRecBestRes");
// mapHisto_["hRecBestRes_Acc"] = new HParticle ("hRecBestRes_Acc");
// mapHisto_["hRecBestResVSMu"] = new HMassVSPart ("hRecBestResVSMu");
//Resolution VS muon kinematic
//----------------------------
mapHisto_["PtResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionGenVSMu");
mapHisto_["PtResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionSimVSMu");
mapHisto_["EtaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionGenVSMu");
mapHisto_["EtaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionSimVSMu");
mapHisto_["ThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionGenVSMu");
mapHisto_["ThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionSimVSMu");
mapHisto_["CotgThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionSimVSMu");
mapHisto_["PhiResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
mapHisto_["PhiResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionSimVSMu");
// Covariances between muons kinematic quantities
// ----------------------------------------------
double ptMax = ptMax_;
// Mass resolution
// ---------------
mapHisto_["MassResolution"] = new HMassResolutionVSPart (outputFile_,"MassResolution");
// mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
// Mass resolution vs (pt, eta) of the muons from MC
massResolutionVsPtEta_ = new HCovarianceVSxy ( "Mass", "Mass", 100, 0., ptMax, 60, -3, 3 );
// Mass resolution vs (pt, eta) of the muons from function
recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
// Muons resolutions from resolution functions
// -------------------------------------------
int totBinsY = 60;
mapHisto_["hFunctionResolMass"] = new HFunctionResolution (outputFile_, "hFunctionResolMass", ptMax, totBinsY);
mapHisto_["hFunctionResolPt"] = new HFunctionResolution (outputFile_, "hFunctionResolPt", ptMax, totBinsY);
mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution (outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
mapHisto_["hFunctionResolPhi"] = new HFunctionResolution (outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
if( readCovariances_ ) {
// Covariances read from file. Used to compare the terms in the expression of mass error
mapHisto_["ReadCovariances"] = new HCovarianceVSParts ( theCovariancesRootFileName_, "Covariance" );
// Variances
mapHisto_["MassResolutionPt"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPt", ptMax);
mapHisto_["MassResolutionCotgTheta"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassCotgTheta", ptMax);
mapHisto_["MassResolutionPhi"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPhi", ptMax);
// Covariances
mapHisto_["MassResolutionPt-CotgTheta"] = new HFunctionResolution(outputFile_,"functionResolMassPt-CotgTheta", ptMax, totBinsY);
mapHisto_["MassResolutionPt-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassPt-Phi", ptMax, totBinsY);
mapHisto_["MassResolutionCotgTheta-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta-Phi", ptMax, totBinsY);
mapHisto_["MassResolutionPt1-Pt2"] = new HFunctionResolution(outputFile_,"functionResolMassPt1-Pt2", ptMax, totBinsY);
mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
mapHisto_["MassResolutionPhi1-Phi2"] = new HFunctionResolution(outputFile_,"functionResolMassPhi1-Phi2", ptMax, totBinsY);
mapHisto_["MassResolutionPt12-CotgTheta21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
mapHisto_["MassResolutionPt12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-Phi21", ptMax, totBinsY);
mapHisto_["MassResolutionCotgTheta12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
mapHisto_["sigmaPtFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPtFromVariance", ptMax, totBinsY);
mapHisto_["sigmaCotgThetaFromVariance"] = new HFunctionResolution(outputFile_,"sigmaCotgThetaFromVariance", ptMax, totBinsY);
mapHisto_["sigmaPhiFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPhiFromVariance", ptMax, totBinsY);
// Derivatives
mapHisto_["DerivativePt"] = new HFunctionResolution(outputFile_,"derivativePt", ptMax);
mapHisto_["DerivativeCotgTheta"] = new HFunctionResolution(outputFile_,"derivativeCotgTheta", ptMax);
mapHisto_["DerivativePhi"] = new HFunctionResolution(outputFile_,"derivativePhi", ptMax);
// Pt term from function
mapHisto_["MassResolutionPtFromFunction"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPtFromFunction", ptMax);
mapHisto_["FullMassResolution"] = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
mapHisto_["MassRes1"] = new HFunctionResolution(outputFile_, "massRes1", ptMax);
mapHisto_["MassRes2"] = new HFunctionResolution(outputFile_, "massRes2", ptMax);
mapHisto_["MassRes3"] = new HFunctionResolution(outputFile_, "massRes3", ptMax);
mapHisto_["MassRes4"] = new HFunctionResolution(outputFile_, "massRes4", ptMax);
mapHisto_["MassRes5"] = new HFunctionResolution(outputFile_, "massRes5", ptMax);
mapHisto_["MassRes6"] = new HFunctionResolution(outputFile_, "massRes6", ptMax);
mapHisto_["MassResPtAndPt12"] = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
}
else {
mapHisto_["Covariances"] = new HCovarianceVSParts ( outputFile_, "Covariance", ptMax );
}
}
| std::vector<reco::LeafCandidate> ResolutionAnalyzer::fillMuonCollection | ( | const std::vector< T > & | tracks | ) | [inline, private] |
Definition at line 76 of file ResolutionAnalyzer.h.
References gather_cfg::cout, debug_, MuScleFitUtils::goodmuon, MuScleFitUtils::mMu2, RPCpg::mu, metsig::muon, patZpeak::muons, and mathSSE::sqrt().
{
std::vector<reco::LeafCandidate> muons;
typename std::vector<T>::const_iterator track;
for (track = tracks.begin(); track != tracks.end(); ++track){
reco::Particle::LorentzVector mu(track->px(),track->py(),track->pz(),
sqrt(track->p()*track->p() + MuScleFitUtils::mMu2));
MuScleFitUtils::goodmuon++;
if (debug_>0) std::cout <<std::setprecision(9)<< "Muon #" << MuScleFitUtils::goodmuon
<< ": initial value Pt = " << mu.Pt() << std::endl;
reco::LeafCandidate muon(track->charge(),mu);
// Store muon
// ----------
muons.push_back( muon );
}
return muons;
}
| void ResolutionAnalyzer::writeHistoMap | ( | ) | [private] |
Writes the histograms in the map.
Definition at line 539 of file ResolutionAnalyzer.cc.
References deltaPtOverPt_, deltaPtOverPtForEta12_, interpolateCardsSimple::histo, mapHisto_, massResolutionVsPtEta_, outputFile_, recoPtVsgenPt_, recoPtVsgenPtEta12_, and HCovarianceVSxy::Write().
Referenced by ~ResolutionAnalyzer().
{
for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin();
histo!=mapHisto_.end(); histo++) {
(*histo).second->Write();
}
outputFile_->cd();
massResolutionVsPtEta_->Write();
recoPtVsgenPt_->Write();
recoPtVsgenPtEta12_->Write();
deltaPtOverPt_->Write();
deltaPtOverPtForEta12_->Write();
}
bool ResolutionAnalyzer::debug_ [private] |
Definition at line 109 of file ResolutionAnalyzer.h.
Referenced by checkDeltaR(), and fillMuonCollection().
TH1D* ResolutionAnalyzer::deltaPtOverPt_ [private] |
Definition at line 125 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
TH1D* ResolutionAnalyzer::deltaPtOverPtForEta12_ [private] |
Definition at line 126 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
int ResolutionAnalyzer::eventCounter_ [private] |
Definition at line 113 of file ResolutionAnalyzer.h.
Referenced by analyze(), ResolutionAnalyzer(), and ~ResolutionAnalyzer().
std::map<std::string, Histograms*> ResolutionAnalyzer::mapHisto_ [private] |
Definition at line 110 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
Definition at line 122 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
int32_t ResolutionAnalyzer::maxEvents_ [private] |
Definition at line 118 of file ResolutionAnalyzer.h.
Referenced by analyze().
TFile* ResolutionAnalyzer::outputFile_ [private] |
Definition at line 111 of file ResolutionAnalyzer.h.
Referenced by fillHistoMap(), ResolutionAnalyzer(), writeHistoMap(), and ~ResolutionAnalyzer().
double ResolutionAnalyzer::ptMax_ [private] |
Definition at line 120 of file ResolutionAnalyzer.h.
Referenced by fillHistoMap().
bool ResolutionAnalyzer::readCovariances_ [private] |
Definition at line 115 of file ResolutionAnalyzer.h.
Referenced by analyze(), and fillHistoMap().
TH2D* ResolutionAnalyzer::recoPtVsgenPt_ [private] |
Definition at line 123 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
TH2D* ResolutionAnalyzer::recoPtVsgenPtEta12_ [private] |
Definition at line 124 of file ResolutionAnalyzer.h.
Referenced by analyze(), fillHistoMap(), and writeHistoMap().
bool ResolutionAnalyzer::resonance_ [private] |
Definition at line 114 of file ResolutionAnalyzer.h.
Referenced by analyze().
std::string ResolutionAnalyzer::theCovariancesRootFileName_ [private] |
Definition at line 108 of file ResolutionAnalyzer.h.
Referenced by fillHistoMap().
Definition at line 104 of file ResolutionAnalyzer.h.
int ResolutionAnalyzer::theMuonType_ [private] |
Definition at line 106 of file ResolutionAnalyzer.h.
std::string ResolutionAnalyzer::theRootFileName_ [private] |
Definition at line 107 of file ResolutionAnalyzer.h.
Referenced by ResolutionAnalyzer().
TString ResolutionAnalyzer::treeFileName_ [private] |
Definition at line 117 of file ResolutionAnalyzer.h.
Referenced by analyze().
1.7.3