CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ResolutionAnalyzer.cc
Go to the documentation of this file.
1 #ifndef RESOLUTIONANALYZER_CC
2 #define RESOLUTIONANALYZER_CC
3 
4 #include "ResolutionAnalyzer.h"
5 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
6 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
7 
8 //
9 // constants, enums and typedefs
10 //
11 
12 //
13 // static data member definitions
14 //
15 
16 //
17 // constructors and destructor
18 //
19 ResolutionAnalyzer::ResolutionAnalyzer(const edm::ParameterSet& iConfig) :
20  theMuonLabel_( iConfig.getParameter<edm::InputTag>( "MuonLabel" ) ),
21  theMuonType_( iConfig.getParameter<int>( "MuonType" ) ),
22  theRootFileName_( iConfig.getUntrackedParameter<std::string>("OutputFileName") ),
23  theCovariancesRootFileName_( iConfig.getUntrackedParameter<std::string>("InputFileName") ),
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];
38  MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType );
39  // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType];
40  MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService( 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 }
52 
53 
54 ResolutionAnalyzer::~ResolutionAnalyzer()
55 {
56  outputFile_->cd();
57  writeHistoMap();
58  outputFile_->Close();
59  std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
60 }
61 
62 
63 //
64 // member functions
65 //
66 
67 // ------------ method called to for each event ------------
68 void ResolutionAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
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 }
419 
420 void ResolutionAnalyzer::fillHistoMap() {
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 }
538 
539 void ResolutionAnalyzer::writeHistoMap() {
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();
545  massResolutionVsPtEta_->Write();
546  recoPtVsgenPt_->Write();
547  recoPtVsgenPtEta12_->Write();
548  deltaPtOverPt_->Write();
549  deltaPtOverPtForEta12_->Write();
550 }
551 
552 bool ResolutionAnalyzer::checkDeltaR(const reco::Particle::LorentzVector & genMu, const reco::Particle::LorentzVector & recMu){
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 }
564 
565 //define this as a plug-in
566 DEFINE_FWK_MODULE(ResolutionAnalyzer);
567 
568 #endif // RESOLUTIONANALYZER_CC
MuScleFitUtils::deltaPhiNoFabs
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...
Definition: MuScleFitUtils.h:101
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
ResolutionAnalyzer::mapHisto_
std::map< std::string, Histograms * > mapHisto_
Definition: ResolutionAnalyzer.h:109
ResolutionAnalyzer::theRootFileName_
std::string theRootFileName_
Definition: ResolutionAnalyzer.h:106
MuScleFitUtils::parResol
static std::vector< double > parResol
Definition: MuScleFitUtils.h:175
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
HCovarianceVSxy::Fill
virtual void Fill(const double &x, const double &y, const double &a, const double &b)
Definition: Histograms.h:1893
ResolutionAnalyzer::massResolutionVsPtEta_
HCovarianceVSxy * massResolutionVsPtEta_
Definition: ResolutionAnalyzer.h:121
ResolutionAnalyzer::fillHistoMap
void fillHistoMap()
Used to fill the map with the histograms needed.
Definition: ResolutionAnalyzer.cc:420
iEvent
int iEvent
Definition: GenABIO.cc:230
resolutionFunctionVecService
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:156
MuScleFitUtils::massResolution
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
ResolutionAnalyzer::writeHistoMap
void writeHistoMap()
Writes the histograms in the map.
Definition: ResolutionAnalyzer.cc:539
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
MuonPairVector
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
Definition: RootTreeHandler.h:13
ResolutionAnalyzer::theCovariancesRootFileName_
std::string theCovariancesRootFileName_
Definition: ResolutionAnalyzer.h:107
ResolutionAnalyzer::checkDeltaR
bool checkDeltaR(const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recMu)
Returns true if the two particles have DeltaR &lt; cut.
Definition: ResolutionAnalyzer.cc:552
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
ResolutionAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: ResolutionAnalyzer.cc:68
MuScleFitUtils::ResolFitType
static int ResolFitType
Definition: MuScleFitUtils.h:135
MuScleFitUtils::mMu2
static const double mMu2
Definition: MuScleFitUtils.h:122
resolutionFunctionService
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:116
MuScleFitUtils::deltaPhi
static double deltaPhi(const double &phi1, const double &phi2)
Definition: MuScleFitUtils.h:93
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
MuScleFitUtils::resfind
static std::vector< int > resfind
Definition: MuScleFitUtils.h:193
MuScleFitUtils::resolutionFunctionForVec
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
Definition: MuScleFitUtils.h:137
reco::Particle::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:27
HResolutionVSPart
A set of histograms for resolution.
Definition: Histograms.h:1139
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
MuScleFitUtils::resolutionFunction
static resolutionFunctionBase< double * > * resolutionFunction
Definition: MuScleFitUtils.h:136
ResolutionAnalyzer::ResolutionAnalyzer
ResolutionAnalyzer(const edm::ParameterSet &)
Definition: ResolutionAnalyzer.cc:19
Generated for CMSSW Reference Manual by doxygen 1.8.5