CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
ResolutionAnalyzer Class Reference

#include <MuonAnalysis/MomentumScaleCalibration/plugins/ResolutionAnalyzer.cc>

Inheritance diagram for ResolutionAnalyzer:
edm::EDAnalyzer edm::EDConsumerBase

Public Member Functions

 ResolutionAnalyzer (const edm::ParameterSet &)
 
 ~ResolutionAnalyzer () override
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
SerialTaskQueueglobalLuminosityBlocksQueue ()
 
SerialTaskQueueglobalRunsQueue ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
 ~EDAnalyzer () override
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
void convertCurrentProcessAlias (std::string const &processName)
 Convert "@currentProcess" in InputTag process names to the actual current process name. More...
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector< ProductResolverIndexAndSkipBit > const & itemsToGetFrom (BranchType iType) const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBaseoperator= (EDConsumerBase &&)=default
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
ProductResolverIndexAndSkipBit uncheckedIndexFrom (EDGetToken) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Member Functions

void analyze (const edm::Event &, const edm::EventSetup &) override
 
bool checkDeltaR (const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recMu)
 Returns true if the two particles have DeltaR < cut. More...
 
void endJob () override
 
void fillHistoMap ()
 Used to fill the map with the histograms needed. More...
 
template<typename T >
std::vector< reco::LeafCandidatefillMuonCollection (const std::vector< T > &tracks)
 
void writeHistoMap ()
 Writes the histograms in the map. More...
 

Private Attributes

bool debug_
 
TH1D * deltaPtOverPt_
 
TH1D * deltaPtOverPtForEta12_
 
int eventCounter_
 
std::map< std::string, Histograms * > mapHisto_
 
HCovarianceVSxymassResolutionVsPtEta_
 
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_
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
static bool wantsGlobalLuminosityBlocks ()
 
static bool wantsGlobalRuns ()
 
static bool wantsStreamLuminosityBlocks ()
 
static bool wantsStreamRuns ()
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes ()
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes (ESInputTag const &tag)
 
template<typename ESProduct , Transition Tr = Transition::Event>
auto esConsumes (eventsetup::EventSetupRecordKey const &, ESInputTag const &tag)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 65 of file ResolutionAnalyzer.h.

Constructor & Destructor Documentation

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_.

19  :
20  theMuonLabel_( iConfig.getParameter<edm::InputTag>( "MuonLabel" ) ),
21  theMuonType_( iConfig.getParameter<int>( "MuonType" ) ),
22  theRootFileName_( iConfig.getUntrackedParameter<std::string>("OutputFileName") ),
24  debug_( iConfig.getUntrackedParameter<bool>( "Debug" ) ),
25  resonance_( iConfig.getParameter<bool>("Resonance") ),
26  readCovariances_( iConfig.getParameter<bool>( "ReadCovariances" ) ),
27  treeFileName_( iConfig.getParameter<std::string>("InputTreeName") ),
28  maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
29  ptMax_( iConfig.getParameter<double>("PtMax") )
30 {
31  //now do what ever initialization is needed
32 
33  // Initial parameters values
34  // -------------------------
35  int resolFitType = iConfig.getParameter<int>("ResolFitType");
36  MuScleFitUtils::ResolFitType = resolFitType;
37  // MuScleFitUtils::resolutionFunction = resolutionFunctionArray[resolFitType];
39  // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType];
41 
42  MuScleFitUtils::parResol = iConfig.getParameter<std::vector<double> >("parResol");
43 
44  MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("ResFind");
45 
46  outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
47  outputFile_->cd();
48  fillHistoMap();
49 
50  eventCounter_ = 0;
51 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
std::string theRootFileName_
static std::vector< double > parResol
void fillHistoMap()
Used to fill the map with the histograms needed.
resolutionFunctionBase< std::vector< double > > * resolutionFunctionVecService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier when receiving a std::...
Definition: Functions.cc:50
std::string theCovariancesRootFileName_
static int ResolFitType
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:38
static std::vector< int > resfind
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
static resolutionFunctionBase< double * > * resolutionFunction
edm::InputTag theMuonLabel_
ResolutionAnalyzer::~ResolutionAnalyzer ( )
override

Definition at line 54 of file ResolutionAnalyzer.cc.

References gather_cfg::cout, eventCounter_, outputFile_, and writeHistoMap().

55 {
56  outputFile_->cd();
57  writeHistoMap();
58  outputFile_->Close();
59  std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
60 }
void writeHistoMap()
Writes the histograms in the map.

Member Function Documentation

void ResolutionAnalyzer::analyze ( const edm::Event iEvent,
const edm::EventSetup iSetup 
)
overrideprivate

Definition at line 68 of file ResolutionAnalyzer.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, checkDeltaR(), funct::cos(), gather_cfg::cout, MuScleFitUtils::deltaPhiNoFabs(), deltaPtOverPt_, deltaPtOverPtForEta12_, reco::tau::disc::Eta(), eventCounter_, JetChargeProducer_cfi::exp, HCovarianceVSxy::Fill(), mps_fire::i, mapHisto_, ResonanceBuilder::mass, MuScleFitUtils::massResolution(), massResolutionVsPtEta_, maxEvents_, MuScleFitUtils::mMu2, MuScleFitUtils::parResol, funct::pow(), reco::tau::disc::Pt(), hiDetachedQuadStep_cff::pt1, hiDetachedQuadStep_cff::pt2, readCovariances_, RootTreeHandler::readTree(), recoPtVsgenPt_, recoPtVsgenPtEta12_, MuScleFitUtils::resolutionFunctionForVec, resonance_, funct::sin(), mathSSE::sqrt(), and treeFileName_.

68  {
69  using namespace edm;
70 
71  std::cout << "starting" << std::endl;
72 
73  lorentzVector nullLorentzVector(0, 0, 0, 0);
74 
75  RootTreeHandler rootTreeHandler;
76  typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
77  MuonPairVector savedPairVector;
78  MuonPairVector genPairVector;
79 
80  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
81  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &evtRun, &genPairVector);
82  MuonPairVector::iterator savedPair = savedPairVector.begin();
83  MuonPairVector::iterator genPair = genPairVector.begin();
84  std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
85  for( ; savedPair != savedPairVector.end(); ++savedPair, ++genPair ) {
86 
87  ++eventCounter_;
88 
89  if( (eventCounter_ % 10000) == 0 ) {
90  std::cout << "event = " << eventCounter_ << std::endl;
91  }
92 
93  lorentzVector recMu1( savedPair->first );
94  lorentzVector recMu2( savedPair->second );
95 
96  if ( resonance_ ) {
97 
98  // Histograms with genParticles characteristics
99  // --------------------------------------------
100 
101  reco::Particle::LorentzVector genMother( genPair->first + genPair->second );
102 
103  mapHisto_["GenMother"]->Fill( genMother );
104  mapHisto_["DeltaGenMotherMuons"]->Fill( genPair->first, genPair->second );
105  mapHisto_["GenMotherMuons"]->Fill( genPair->first );
106  mapHisto_["GenMotherMuons"]->Fill( genPair->second );
107 
108  // Match the reco muons with the gen and sim tracks
109  // ------------------------------------------------
110  if(checkDeltaR(genPair->first,recMu1)){
111  mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Pt()+recMu1.Pt())/genPair->first.Pt(),-1);
112  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Theta()+recMu1.Theta()),-1);
113  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu1,(-cos(genPair->first.Theta())/sin(genPair->first.Theta())
114  +cos(recMu1.Theta())/sin(recMu1.Theta())),-1);
115  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Eta()+recMu1.Eta()),-1);
116  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
117  mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()),-1);
118  recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
119  deltaPtOverPt_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
120  if( fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2 ) {
121  recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
122  deltaPtOverPtForEta12_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
123  }
124  }
125  if(checkDeltaR(genPair->second,recMu2)){
126  mapHisto_["PtResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Pt()+recMu2.Pt())/genPair->second.Pt(),+1);
127  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Theta()+recMu2.Theta()),+1);
128  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu2,(-cos(genPair->second.Theta())/sin(genPair->second.Theta())
129  +cos(recMu2.Theta())/sin(recMu2.Theta())),+1);
130  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Eta()+recMu2.Eta()),+1);
131  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
132  mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()),+1);
133  recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
134  deltaPtOverPt_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
135  if( fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2 ) {
136  recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
137  deltaPtOverPtForEta12_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
138  }
139  }
140 
141  // Fill the mass resolution histograms
142  // -----------------------------------
143  // check if the recoMuons match the genMuons
144  // if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
145  if( genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
146  checkDeltaR(genPair->first,recMu1) && checkDeltaR(genPair->second,recMu2) ) {
147 
148  double recoMass = (recMu1+recMu2).mass();
149  double genMass = (genPair->first + genPair->second).mass();
150  // first is always mu-, second is always mu+
151  mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
152 
153  // Fill the reconstructed resonance
154  reco::Particle::LorentzVector recoResonance( recMu1+recMu2 );
155  mapHisto_["RecoResonance"]->Fill( recoResonance );
156  mapHisto_["DeltaRecoResonanceMuons"]->Fill( recMu1, recMu2 );
157  mapHisto_["RecoResonanceMuons"]->Fill( recMu1 );
158  mapHisto_["RecoResonanceMuons"]->Fill( recMu2 );
159 
160  // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
161  if( genMass != 0 ) {
162  // double diffMass = (recoMass - genMass)/genMass;
163  double diffMass = recoMass - genMass;
164  // Fill if for both muons
165  double pt1 = recMu1.pt();
166  double eta1 = recMu1.eta();
167  double pt2 = recMu2.pt();
168  double eta2 = recMu2.eta();
169  // This is to avoid nan
170  if( diffMass == diffMass ) {
171  massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
172  massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
173  }
174  else {
175  std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
176  }
177  // Fill with mass resolution from resolution function
178  double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
179  // The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
180  // We alos take the squared value, since var = sigma^2.
181  mapHisto_["hFunctionResolMass"]->Fill( recMu1, std::pow(massRes,2), -1 );
182  mapHisto_["hFunctionResolMass"]->Fill( recMu2, std::pow(massRes,2), +1 );
183  }
184 
185  // Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
186  mapHisto_["hFunctionResolPt"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
187  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
188  mapHisto_["hFunctionResolPhi"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
189  mapHisto_["hFunctionResolPt"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
190  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
191  mapHisto_["hFunctionResolPhi"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
192 
193  if( readCovariances_ ) {
194  // Compute mass error terms
195  // ------------------------
196  double mass = (recMu1+recMu2).mass();
197  double pt1 = recMu1.Pt();
198  double phi1 = recMu1.Phi();
199  double eta1 = recMu1.Eta();
200  double theta1 = 2*atan(exp(-eta1));
201  double pt2 = recMu2.Pt();
202  double phi2 = recMu2.Phi();
203  double eta2 = recMu2.Eta();
204  double theta2 = 2*atan(exp(-eta2));
205  // Derivatives
206  double mMu2 = MuScleFitUtils::mMu2;
207  double dmdpt1 = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))-
208  pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
209  double dmdpt2 = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))-
210  pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
211  double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
212  double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
213  double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
214  sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) -
215  pt1*pt2*cos(theta2)/sin(theta2))/mass;
216  double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
217  sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) -
218  pt2*pt1*cos(theta1)/sin(theta1))/mass;
219 
220  // Multiplied by the pt here
221  // -------------------------
222  double dmdpt[2] = {dmdpt1*recMu1.Pt(), dmdpt2*recMu2.Pt()};
223  double dmdphi[2] = {dmdphi1, dmdphi2};
224  double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
225 
226  // Evaluate the single terms in the mass error expression
227 
228  reco::Particle::LorentzVector * recMu[2] = { &recMu1, &recMu2 };
229  int charge[2] = { -1, +1 };
230 
231  double fullMassRes = 0.;
232  double massRes1 = 0.;
233  double massRes2 = 0.;
234  double massRes3 = 0.;
235  double massRes4 = 0.;
236  double massRes5 = 0.;
237  double massRes6 = 0.;
238  double massResPtAndPt12 = 0.;
239 
240  for( int i=0; i<2; ++i ) {
241 
242  double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
243  double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
244  double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
245  double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
246  double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
247  double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
248 
249  double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
250  double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
251  double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
252  double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
253  double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
254  double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
255 
256  // ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
257  mapHisto_["MassResolutionPt"]->Fill( *(recMu[i]), ptVariance*std::pow(dmdpt[i],2), charge[i] );
258  mapHisto_["MassResolutionCotgTheta"]->Fill( *(recMu[i]), cotgThetaVariance*std::pow(dmdcotgth[i],2), charge[i] );
259  mapHisto_["MassResolutionPhi"]->Fill( *(recMu[i]), phiVariance*std::pow(dmdphi[i],2), charge[i] );
260  mapHisto_["MassResolutionPt-CotgTheta"]->Fill( *(recMu[i]), 2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i], charge[i] );
261  mapHisto_["MassResolutionPt-Phi"]->Fill( *(recMu[i]), 2*pt_phi*dmdpt[i]*dmdphi[i], charge[i] );
262  mapHisto_["MassResolutionCotgTheta-Phi"]->Fill( *(recMu[i]), 2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i], charge[i] );
263 
264  mapHisto_["MassResolutionPt1-Pt2"]->Fill( *(recMu[i]), pt1_pt2*dmdpt[0]*dmdpt[1], charge[i] );
265  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill( *(recMu[i]), cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1], charge[i] );
266  mapHisto_["MassResolutionPhi1-Phi2"]->Fill( *(recMu[i]), phi1_phi2*dmdphi[0]*dmdphi[1], charge[i] );
267  // This must be filled for both configurations: 12 and 21
268  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1], charge[i] );
269  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0], charge[i] );
270  mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[0]*dmdphi[1], charge[i] );
271  mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[1]*dmdphi[0], charge[i] );
272  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1], charge[i] );
273  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0], charge[i] );
274 
275  // Sigmas for comparison
276  mapHisto_["sigmaPtFromVariance"]->Fill( *(recMu[i]), sqrt(ptVariance), charge[i] );
277  mapHisto_["sigmaCotgThetaFromVariance"]->Fill( *(recMu[i]), sqrt(cotgThetaVariance), charge[i] );
278  mapHisto_["sigmaPhiFromVariance"]->Fill( *(recMu[i]), sqrt(phiVariance), charge[i] );
279 
280  // Pt term from function
281  mapHisto_["MassResolutionPtFromFunction"]->Fill( *(recMu[i]), ( MuScleFitUtils::resolutionFunctionForVec->sigmaPt((recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol) )*std::pow(dmdpt[i],2), charge[i] );
282 
283  fullMassRes +=
284  ptVariance*std::pow(dmdpt[i],2) +
285  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
286  phiVariance*std::pow(dmdphi[i],2) +
287 
288  // These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
289  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
290  2*pt_phi*dmdpt[i]*dmdphi[i] +
291  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
292 
293  pt1_pt2*dmdpt[0]*dmdpt[1] +
294  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
295  phi1_phi2*dmdphi[0]*dmdphi[1] +
296 
297  // These are filled twice, because of the two combinations
298  pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
299  pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
300  pt12_phi21*dmdpt[0]*dmdphi[1] +
301  pt12_phi21*dmdpt[1]*dmdphi[0] +
302  cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
303  cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
304 
305  massRes1 += ptVariance*std::pow(dmdpt[i],2);
306  massRes2 += ptVariance*std::pow(dmdpt[i],2) +
307  cotgThetaVariance*std::pow(dmdcotgth[i],2);
308  massRes3 += ptVariance*std::pow(dmdpt[i],2) +
309  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
310  phiVariance*std::pow(dmdphi[i],2);
311  massRes4 += ptVariance*std::pow(dmdpt[i],2) +
312  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
313  phiVariance*std::pow(dmdphi[i],2) +
314  pt1_pt2*dmdpt[0]*dmdpt[1] +
315  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
316  2*pt_phi*dmdpt[i]*dmdphi[i] +
317  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i];
318  massRes5 += ptVariance*std::pow(dmdpt[i],2) +
319  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
320  phiVariance*std::pow(dmdphi[i],2) +
321  pt1_pt2*dmdpt[0]*dmdpt[1] +
322  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
323  2*pt_phi*dmdpt[i]*dmdphi[i] +
324  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
325  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
326  phi1_phi2*dmdphi[0]*dmdphi[1];
327  massRes6 += ptVariance*std::pow(dmdpt[i],2) +
328  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
329  phiVariance*std::pow(dmdphi[i],2) +
330  pt1_pt2*dmdpt[0]*dmdpt[1] +
331  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
332  2*pt_phi*dmdpt[i]*dmdphi[i] +
333  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
334  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
335  phi1_phi2*dmdphi[0]*dmdphi[1] +
336  pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
337  pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
338  pt12_phi21*dmdpt[0]*dmdphi[1] +
339  pt12_phi21*dmdpt[1]*dmdphi[0] +
340  cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
341  cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
342 
343  massResPtAndPt12 += ptVariance*std::pow(dmdpt[i],2) + pt1_pt2*dmdpt[0]*dmdpt[1];
344 
345  // Derivatives
346  mapHisto_["DerivativePt"]->Fill( *(recMu[i]), dmdpt[i], charge[i] );
347  mapHisto_["DerivativeCotgTheta"]->Fill( *(recMu[i]), dmdcotgth[i], charge[i] );
348  mapHisto_["DerivativePhi"]->Fill( *(recMu[i]), dmdphi[i], charge[i] );
349  }
350  // Fill the complete resolution function with covariance terms
351  mapHisto_["FullMassResolution"]->Fill( *(recMu[0]), fullMassRes, charge[0]);
352  mapHisto_["FullMassResolution"]->Fill( *(recMu[1]), fullMassRes, charge[1]);
353 
354  mapHisto_["MassRes1"]->Fill( *(recMu[0]), massRes1, charge[0] );
355  mapHisto_["MassRes1"]->Fill( *(recMu[1]), massRes1, charge[1] );
356  mapHisto_["MassRes2"]->Fill( *(recMu[0]), massRes2, charge[0] );
357  mapHisto_["MassRes2"]->Fill( *(recMu[1]), massRes2, charge[1] );
358  mapHisto_["MassRes3"]->Fill( *(recMu[0]), massRes3, charge[0] );
359  mapHisto_["MassRes3"]->Fill( *(recMu[1]), massRes3, charge[1] );
360  mapHisto_["MassRes4"]->Fill( *(recMu[0]), massRes4, charge[0] );
361  mapHisto_["MassRes4"]->Fill( *(recMu[1]), massRes4, charge[1] );
362  mapHisto_["MassRes5"]->Fill( *(recMu[0]), massRes5, charge[0] );
363  mapHisto_["MassRes5"]->Fill( *(recMu[1]), massRes5, charge[1] );
364  mapHisto_["MassRes6"]->Fill( *(recMu[0]), massRes6, charge[0] );
365  mapHisto_["MassRes6"]->Fill( *(recMu[1]), massRes6, charge[1] );
366  mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[0]), massResPtAndPt12, charge[0] );
367  mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[1]), massResPtAndPt12, charge[1] );
368  }
369  else {
370  // Fill the covariances histograms
371  mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
372  }
373  }
374  } // end if resonance
375  }
376 // else {
377 //
378 // // Loop on the recMuons
379 // std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
380 // for ( ; recMuon!=muons.end(); ++recMuon ) {
381 // int charge = recMuon->charge();
382 //
383 // lorentzVector recMu(recMuon->p4());
384 //
385 // // Find the matching MC muon
386 // const HepMC::GenEvent* Evt = evtMC->GetEvent();
387 // //Loop on generated particles
388 // std::map<double, lorentzVector> genAssocMap;
389 // HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
390 // for( ; part!=Evt->particles_end(); ++part ) {
391 // if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
392 // lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
393 // (*part)->momentum().pz(),(*part)->momentum().e()));
394 //
395 // double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
396 // ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
397 //
398 // // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
399 // // Only in this case we consider it matching.
400 // if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
401 // }
402 // }
403 // // Take the closest in deltaR
404 // lorentzVector genMu(genAssocMap.begin()->second);
405 //
406 // // Histograms with genParticles characteristics
407 // // --------------------------------------------
408 //
409 // if(checkDeltaR(genMu,recMu)){
410 // mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
411 // mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
412 // mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
413 // +cos(recMu.Theta())/sin(recMu.Theta())),charge);
414 // mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
415 // mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
416 // }
417 // }
418 }
static double deltaPhiNoFabs(const double &phi1, const double &phi2)
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance c...
std::map< std::string, Histograms * > mapHisto_
static std::vector< double > parResol
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void Fill(const double &x, const double &y, const double &a, const double &b) override
Definition: Histograms.h:1893
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, std::vector< std::pair< unsigned int, unsigned long long > > *evtRun, MuonPairVector *genPair=nullptr)
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
HCovarianceVSxy * massResolutionVsPtEta_
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
T sqrt(T t)
Definition: SSEVec.h:18
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
bool checkDeltaR(const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recMu)
Returns true if the two particles have DeltaR < cut.
static const double mMu2
HLT enums.
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
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_, DEFINE_FWK_MODULE, MuScleFitUtils::deltaPhi(), deltaR(), and mathSSE::sqrt().

Referenced by analyze(), and fillMuonCollection().

552  {
553  //first is always mu-, second is always mu+
554  double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) +
555  ((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
556  if(deltaR<0.01)
557  return true;
558  else if (debug_>0)
559  std::cout<<"Reco muon "<<recMu<<" with eta "<<recMu.Eta()<<" and phi "<<recMu.Phi()<<std::endl
560  <<" DOES NOT MATCH with generated muon from resonance: "<<std::endl
561  <<genMu<<" with eta "<<genMu.Eta()<<" and phi "<<genMu.Phi()<<std::endl;
562  return false;
563 }
T sqrt(T t)
Definition: SSEVec.h:18
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
static double deltaPhi(const double &phi1, const double &phi2)
void ResolutionAnalyzer::endJob ( void  )
inlineoverrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 72 of file ResolutionAnalyzer.h.

72 {};
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_, PFRecoTauDiscriminationByHPSSelection_cfi::maxMass, PFRecoTauDiscriminationByHPSSelection_cfi::minMass, outputFile_, AlignmentTrackSelector_cfi::ptMax, ptMax_, readCovariances_, recoPtVsgenPt_, recoPtVsgenPtEta12_, MuScleFitUtils::resfind, and theCovariancesRootFileName_.

Referenced by fillMuonCollection(), and ResolutionAnalyzer().

420  {
421 
422  outputFile_->cd();
423 
424  // Resonances
425  // If no Z is required, use a smaller mass range.
426  double minMass = 0.;
427  double maxMass = 200.;
428  if( MuScleFitUtils::resfind[0] != 1 ) {
429  maxMass = 30.;
430  }
431  mapHisto_["GenMother"] = new HParticle(outputFile_, "GenMother", minMass, maxMass);
432  mapHisto_["SimResonance"] = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
433  mapHisto_["RecoResonance"] = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
434 
435  // Resonance muons
436  mapHisto_["GenMotherMuons"] = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
437  mapHisto_["SimResonanceMuons"] = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
438  mapHisto_["RecoResonanceMuons"] = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
439 
440  // Deltas between resonance muons
441  mapHisto_["DeltaGenMotherMuons"] = new HDelta (outputFile_, "DeltaGenMotherMuons");
442  mapHisto_["DeltaSimResonanceMuons"] = new HDelta (outputFile_, "DeltaSimResonanceMuons");
443  mapHisto_["DeltaRecoResonanceMuons"] = new HDelta (outputFile_, "DeltaRecoResonanceMuons");
444 
445  // //Reconstructed muon kinematics
446  // //-----------------------------
447  // mapHisto_["hRecBestMu"] = new HParticle ("hRecBestMu");
448  // mapHisto_["hRecBestMu_Acc"] = new HParticle ("hRecBestMu_Acc");
449  // mapHisto_["hDeltaRecBestMu"] = new HDelta ("hDeltaRecBestMu");
450 
451  // mapHisto_["hRecBestRes"] = new HParticle ("hRecBestRes");
452  // mapHisto_["hRecBestRes_Acc"] = new HParticle ("hRecBestRes_Acc");
453  // mapHisto_["hRecBestResVSMu"] = new HMassVSPart ("hRecBestResVSMu");
454 
455  //Resolution VS muon kinematic
456  //----------------------------
457  mapHisto_["PtResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionGenVSMu");
458  mapHisto_["PtResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionSimVSMu");
459  mapHisto_["EtaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionGenVSMu");
460  mapHisto_["EtaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionSimVSMu");
461  mapHisto_["ThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionGenVSMu");
462  mapHisto_["ThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionSimVSMu");
463  mapHisto_["CotgThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
464  mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionSimVSMu");
465  mapHisto_["PhiResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
466  mapHisto_["PhiResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionSimVSMu");
467 
468  // Covariances between muons kinematic quantities
469  // ----------------------------------------------
470  double ptMax = ptMax_;
471 
472  // Mass resolution
473  // ---------------
474  mapHisto_["MassResolution"] = new HMassResolutionVSPart (outputFile_,"MassResolution");
475 
476  // mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
477 
478  // Mass resolution vs (pt, eta) of the muons from MC
479  massResolutionVsPtEta_ = new HCovarianceVSxy ( "Mass", "Mass", 100, 0., ptMax, 60, -3, 3 );
480  // Mass resolution vs (pt, eta) of the muons from function
481  recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
482  recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
483  deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
484  deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
485 
486  // Muons resolutions from resolution functions
487  // -------------------------------------------
488  int totBinsY = 60;
489  mapHisto_["hFunctionResolMass"] = new HFunctionResolution (outputFile_, "hFunctionResolMass", ptMax, totBinsY);
490  mapHisto_["hFunctionResolPt"] = new HFunctionResolution (outputFile_, "hFunctionResolPt", ptMax, totBinsY);
491  mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution (outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
492  mapHisto_["hFunctionResolPhi"] = new HFunctionResolution (outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
493 
494  if( readCovariances_ ) {
495  // Covariances read from file. Used to compare the terms in the expression of mass error
496  mapHisto_["ReadCovariances"] = new HCovarianceVSParts ( theCovariancesRootFileName_, "Covariance" );
497 
498  // Variances
499  mapHisto_["MassResolutionPt"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPt", ptMax);
500  mapHisto_["MassResolutionCotgTheta"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassCotgTheta", ptMax);
501  mapHisto_["MassResolutionPhi"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPhi", ptMax);
502  // Covariances
503  mapHisto_["MassResolutionPt-CotgTheta"] = new HFunctionResolution(outputFile_,"functionResolMassPt-CotgTheta", ptMax, totBinsY);
504  mapHisto_["MassResolutionPt-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassPt-Phi", ptMax, totBinsY);
505  mapHisto_["MassResolutionCotgTheta-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta-Phi", ptMax, totBinsY);
506  mapHisto_["MassResolutionPt1-Pt2"] = new HFunctionResolution(outputFile_,"functionResolMassPt1-Pt2", ptMax, totBinsY);
507  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
508  mapHisto_["MassResolutionPhi1-Phi2"] = new HFunctionResolution(outputFile_,"functionResolMassPhi1-Phi2", ptMax, totBinsY);
509  mapHisto_["MassResolutionPt12-CotgTheta21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
510  mapHisto_["MassResolutionPt12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-Phi21", ptMax, totBinsY);
511  mapHisto_["MassResolutionCotgTheta12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
512 
513  mapHisto_["sigmaPtFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPtFromVariance", ptMax, totBinsY);
514  mapHisto_["sigmaCotgThetaFromVariance"] = new HFunctionResolution(outputFile_,"sigmaCotgThetaFromVariance", ptMax, totBinsY);
515  mapHisto_["sigmaPhiFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPhiFromVariance", ptMax, totBinsY);
516 
517  // Derivatives
518  mapHisto_["DerivativePt"] = new HFunctionResolution(outputFile_,"derivativePt", ptMax);
519  mapHisto_["DerivativeCotgTheta"] = new HFunctionResolution(outputFile_,"derivativeCotgTheta", ptMax);
520  mapHisto_["DerivativePhi"] = new HFunctionResolution(outputFile_,"derivativePhi", ptMax);
521 
522  // Pt term from function
523  mapHisto_["MassResolutionPtFromFunction"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPtFromFunction", ptMax);
524 
525  mapHisto_["FullMassResolution"] = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
526  mapHisto_["MassRes1"] = new HFunctionResolution(outputFile_, "massRes1", ptMax);
527  mapHisto_["MassRes2"] = new HFunctionResolution(outputFile_, "massRes2", ptMax);
528  mapHisto_["MassRes3"] = new HFunctionResolution(outputFile_, "massRes3", ptMax);
529  mapHisto_["MassRes4"] = new HFunctionResolution(outputFile_, "massRes4", ptMax);
530  mapHisto_["MassRes5"] = new HFunctionResolution(outputFile_, "massRes5", ptMax);
531  mapHisto_["MassRes6"] = new HFunctionResolution(outputFile_, "massRes6", ptMax);
532  mapHisto_["MassResPtAndPt12"] = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
533  }
534  else {
535  mapHisto_["Covariances"] = new HCovarianceVSParts ( outputFile_, "Covariance", ptMax );
536  }
537 }
std::map< std::string, Histograms * > mapHisto_
HCovarianceVSxy * massResolutionVsPtEta_
std::string theCovariancesRootFileName_
static std::vector< int > resfind
A set of histograms for resolution.
Definition: Histograms.h:1139
template<typename T >
std::vector<reco::LeafCandidate> ResolutionAnalyzer::fillMuonCollection ( const std::vector< T > &  tracks)
inlineprivate

Definition at line 75 of file ResolutionAnalyzer.h.

References checkDeltaR(), gather_cfg::cout, debug_, fillHistoMap(), MuScleFitUtils::goodmuon, MuScleFitUtils::mMu2, RPCpg::mu, metsig::muon, extraflags_cff::muons, mathSSE::sqrt(), HiIsolationCommonParameters_cff::track, and writeHistoMap().

75  {
76  std::vector<reco::LeafCandidate> muons;
77  typename std::vector<T>::const_iterator track;
78  for (track = tracks.begin(); track != tracks.end(); ++track){
79  reco::Particle::LorentzVector mu(track->px(),track->py(),track->pz(),
80  sqrt(track->p()*track->p() + MuScleFitUtils::mMu2));
82  if (debug_>0) std::cout <<std::setprecision(9)<< "Muon #" << MuScleFitUtils::goodmuon
83  << ": initial value Pt = " << mu.Pt() << std::endl;
84  reco::LeafCandidate muon(track->charge(),mu);
85  // Store muon
86  // ----------
87  muons.push_back( muon );
88  }
89  return muons;
90  }
T sqrt(T t)
Definition: SSEVec.h:18
const int mu
Definition: Constants.h:22
static int goodmuon
static const double mMu2
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
void ResolutionAnalyzer::writeHistoMap ( )
private

Writes the histograms in the map.

Definition at line 539 of file ResolutionAnalyzer.cc.

References deltaPtOverPt_, deltaPtOverPtForEta12_, trackerHits::histo, mapHisto_, massResolutionVsPtEta_, outputFile_, recoPtVsgenPt_, recoPtVsgenPtEta12_, and HCovarianceVSxy::Write().

Referenced by fillMuonCollection(), and ~ResolutionAnalyzer().

539  {
540  for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin();
541  histo!=mapHisto_.end(); histo++) {
542  (*histo).second->Write();
543  }
544  outputFile_->cd();
546  recoPtVsgenPt_->Write();
547  recoPtVsgenPtEta12_->Write();
548  deltaPtOverPt_->Write();
549  deltaPtOverPtForEta12_->Write();
550 }
std::map< std::string, Histograms * > mapHisto_
HCovarianceVSxy * massResolutionVsPtEta_
void Write() override
Definition: Histograms.h:1918

Member Data Documentation

bool ResolutionAnalyzer::debug_
private

Definition at line 108 of file ResolutionAnalyzer.h.

Referenced by checkDeltaR(), and fillMuonCollection().

TH1D* ResolutionAnalyzer::deltaPtOverPt_
private

Definition at line 124 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

TH1D* ResolutionAnalyzer::deltaPtOverPtForEta12_
private

Definition at line 125 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

int ResolutionAnalyzer::eventCounter_
private

Definition at line 112 of file ResolutionAnalyzer.h.

Referenced by analyze(), ResolutionAnalyzer(), and ~ResolutionAnalyzer().

std::map<std::string, Histograms*> ResolutionAnalyzer::mapHisto_
private

Definition at line 109 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

HCovarianceVSxy* ResolutionAnalyzer::massResolutionVsPtEta_
private

Definition at line 121 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

int32_t ResolutionAnalyzer::maxEvents_
private

Definition at line 117 of file ResolutionAnalyzer.h.

Referenced by analyze().

TFile* ResolutionAnalyzer::outputFile_
private
double ResolutionAnalyzer::ptMax_
private

Definition at line 119 of file ResolutionAnalyzer.h.

Referenced by fillHistoMap().

bool ResolutionAnalyzer::readCovariances_
private

Definition at line 114 of file ResolutionAnalyzer.h.

Referenced by analyze(), and fillHistoMap().

TH2D* ResolutionAnalyzer::recoPtVsgenPt_
private

Definition at line 122 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

TH2D* ResolutionAnalyzer::recoPtVsgenPtEta12_
private

Definition at line 123 of file ResolutionAnalyzer.h.

Referenced by analyze(), fillHistoMap(), and writeHistoMap().

bool ResolutionAnalyzer::resonance_
private

Definition at line 113 of file ResolutionAnalyzer.h.

Referenced by analyze().

std::string ResolutionAnalyzer::theCovariancesRootFileName_
private

Definition at line 107 of file ResolutionAnalyzer.h.

Referenced by fillHistoMap().

edm::InputTag ResolutionAnalyzer::theMuonLabel_
private

Definition at line 103 of file ResolutionAnalyzer.h.

int ResolutionAnalyzer::theMuonType_
private

Definition at line 105 of file ResolutionAnalyzer.h.

std::string ResolutionAnalyzer::theRootFileName_
private

Definition at line 106 of file ResolutionAnalyzer.h.

Referenced by ResolutionAnalyzer().

TString ResolutionAnalyzer::treeFileName_
private

Definition at line 116 of file ResolutionAnalyzer.h.

Referenced by analyze().