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  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &genPairVector);
80  MuonPairVector::iterator savedPair = savedPairVector.begin();
81  MuonPairVector::iterator genPair = genPairVector.begin();
82  std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
83  for( ; savedPair != savedPairVector.end(); ++savedPair, ++genPair ) {
84 
85  ++eventCounter_;
86 
87  if( (eventCounter_ % 10000) == 0 ) {
88  std::cout << "event = " << eventCounter_ << std::endl;
89  }
90 
91  lorentzVector recMu1( savedPair->first );
92  lorentzVector recMu2( savedPair->second );
93 
94  if ( resonance_ ) {
95 
96  // Histograms with genParticles characteristics
97  // --------------------------------------------
98 
99  reco::Particle::LorentzVector genMother( genPair->first + genPair->second );
100 
101  mapHisto_["GenMother"]->Fill( genMother );
102  mapHisto_["DeltaGenMotherMuons"]->Fill( genPair->first, genPair->second );
103  mapHisto_["GenMotherMuons"]->Fill( genPair->first );
104  mapHisto_["GenMotherMuons"]->Fill( genPair->second );
105 
106  // Match the reco muons with the gen and sim tracks
107  // ------------------------------------------------
108  if(checkDeltaR(genPair->first,recMu1)){
109  mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Pt()+recMu1.Pt())/genPair->first.Pt(),-1);
110  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Theta()+recMu1.Theta()),-1);
111  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu1,(-cos(genPair->first.Theta())/sin(genPair->first.Theta())
112  +cos(recMu1.Theta())/sin(recMu1.Theta())),-1);
113  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Eta()+recMu1.Eta()),-1);
114  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
115  mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()),-1);
116  recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
117  deltaPtOverPt_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
118  if( fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2 ) {
119  recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
120  deltaPtOverPtForEta12_->Fill( (recMu1.Pt() - genPair->first.Pt())/genPair->first.Pt() );
121  }
122  }
123  if(checkDeltaR(genPair->second,recMu2)){
124  mapHisto_["PtResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Pt()+recMu2.Pt())/genPair->second.Pt(),+1);
125  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Theta()+recMu2.Theta()),+1);
126  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(recMu2,(-cos(genPair->second.Theta())/sin(genPair->second.Theta())
127  +cos(recMu2.Theta())/sin(recMu2.Theta())),+1);
128  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Eta()+recMu2.Eta()),+1);
129  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
130  mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()),+1);
131  recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
132  deltaPtOverPt_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
133  if( fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2 ) {
134  recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
135  deltaPtOverPtForEta12_->Fill( (recMu2.Pt() - genPair->second.Pt())/genPair->second.Pt() );
136  }
137  }
138 
139  // Fill the mass resolution histograms
140  // -----------------------------------
141  // check if the recoMuons match the genMuons
142  // if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
143  if( genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
144  checkDeltaR(genPair->first,recMu1) && checkDeltaR(genPair->second,recMu2) ) {
145 
146  double recoMass = (recMu1+recMu2).mass();
147  double genMass = (genPair->first + genPair->second).mass();
148  // first is always mu-, second is always mu+
149  mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
150 
151  // Fill the reconstructed resonance
152  reco::Particle::LorentzVector recoResonance( recMu1+recMu2 );
153  mapHisto_["RecoResonance"]->Fill( recoResonance );
154  mapHisto_["DeltaRecoResonanceMuons"]->Fill( recMu1, recMu2 );
155  mapHisto_["RecoResonanceMuons"]->Fill( recMu1 );
156  mapHisto_["RecoResonanceMuons"]->Fill( recMu2 );
157 
158  // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
159  if( genMass != 0 ) {
160  // double diffMass = (recoMass - genMass)/genMass;
161  double diffMass = recoMass - genMass;
162  // Fill if for both muons
163  double pt1 = recMu1.pt();
164  double eta1 = recMu1.eta();
165  double pt2 = recMu2.pt();
166  double eta2 = recMu2.eta();
167  // This is to avoid nan
168  if( diffMass == diffMass ) {
169  massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
170  massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
171  }
172  else {
173  std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
174  }
175  // Fill with mass resolution from resolution function
176  double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
177  // The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
178  // We alos take the squared value, since var = sigma^2.
179  mapHisto_["hFunctionResolMass"]->Fill( recMu1, std::pow(massRes,2), -1 );
180  mapHisto_["hFunctionResolMass"]->Fill( recMu2, std::pow(massRes,2), +1 );
181  }
182 
183  // Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
184  mapHisto_["hFunctionResolPt"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
185  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
186  mapHisto_["hFunctionResolPhi"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol), -1 );
187  mapHisto_["hFunctionResolPt"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
188  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
189  mapHisto_["hFunctionResolPhi"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol), +1 );
190 
191  if( readCovariances_ ) {
192  // Compute mass error terms
193  // ------------------------
194  double mass = (recMu1+recMu2).mass();
195  double pt1 = recMu1.Pt();
196  double phi1 = recMu1.Phi();
197  double eta1 = recMu1.Eta();
198  double theta1 = 2*atan(exp(-eta1));
199  double pt2 = recMu2.Pt();
200  double phi2 = recMu2.Phi();
201  double eta2 = recMu2.Eta();
202  double theta2 = 2*atan(exp(-eta2));
203  // Derivatives
204  double mMu2 = MuScleFitUtils::mMu2;
205  double dmdpt1 = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))-
206  pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
207  double dmdpt2 = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))-
208  pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
209  double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
210  double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
211  double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
212  sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) -
213  pt1*pt2*cos(theta2)/sin(theta2))/mass;
214  double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
215  sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) -
216  pt2*pt1*cos(theta1)/sin(theta1))/mass;
217 
218  // Multiplied by the pt here
219  // -------------------------
220  double dmdpt[2] = {dmdpt1*recMu1.Pt(), dmdpt2*recMu2.Pt()};
221  double dmdphi[2] = {dmdphi1, dmdphi2};
222  double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
223 
224  // Evaluate the single terms in the mass error expression
225 
226  reco::Particle::LorentzVector * recMu[2] = { &recMu1, &recMu2 };
227  int charge[2] = { -1, +1 };
228 
229  double fullMassRes = 0.;
230  double massRes1 = 0.;
231  double massRes2 = 0.;
232  double massRes3 = 0.;
233  double massRes4 = 0.;
234  double massRes5 = 0.;
235  double massRes6 = 0.;
236  double massResPtAndPt12 = 0.;
237 
238  for( int i=0; i<2; ++i ) {
239 
240  double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
241  double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
242  double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
243  double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
244  double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
245  double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
246 
247  double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
248  double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
249  double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
250  double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
251  double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
252  double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
253 
254  // ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
255  mapHisto_["MassResolutionPt"]->Fill( *(recMu[i]), ptVariance*std::pow(dmdpt[i],2), charge[i] );
256  mapHisto_["MassResolutionCotgTheta"]->Fill( *(recMu[i]), cotgThetaVariance*std::pow(dmdcotgth[i],2), charge[i] );
257  mapHisto_["MassResolutionPhi"]->Fill( *(recMu[i]), phiVariance*std::pow(dmdphi[i],2), charge[i] );
258  mapHisto_["MassResolutionPt-CotgTheta"]->Fill( *(recMu[i]), 2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i], charge[i] );
259  mapHisto_["MassResolutionPt-Phi"]->Fill( *(recMu[i]), 2*pt_phi*dmdpt[i]*dmdphi[i], charge[i] );
260  mapHisto_["MassResolutionCotgTheta-Phi"]->Fill( *(recMu[i]), 2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i], charge[i] );
261 
262  mapHisto_["MassResolutionPt1-Pt2"]->Fill( *(recMu[i]), pt1_pt2*dmdpt[0]*dmdpt[1], charge[i] );
263  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill( *(recMu[i]), cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1], charge[i] );
264  mapHisto_["MassResolutionPhi1-Phi2"]->Fill( *(recMu[i]), phi1_phi2*dmdphi[0]*dmdphi[1], charge[i] );
265  // This must be filled for both configurations: 12 and 21
266  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1], charge[i] );
267  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill( *(recMu[i]), pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0], charge[i] );
268  mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[0]*dmdphi[1], charge[i] );
269  mapHisto_["MassResolutionPt12-Phi21"]->Fill( *(recMu[i]), pt12_phi21*dmdpt[1]*dmdphi[0], charge[i] );
270  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1], charge[i] );
271  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill( *(recMu[i]), cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0], charge[i] );
272 
273  // Sigmas for comparison
274  mapHisto_["sigmaPtFromVariance"]->Fill( *(recMu[i]), sqrt(ptVariance), charge[i] );
275  mapHisto_["sigmaCotgThetaFromVariance"]->Fill( *(recMu[i]), sqrt(cotgThetaVariance), charge[i] );
276  mapHisto_["sigmaPhiFromVariance"]->Fill( *(recMu[i]), sqrt(phiVariance), charge[i] );
277 
278  // Pt term from function
279  mapHisto_["MassResolutionPtFromFunction"]->Fill( *(recMu[i]), ( MuScleFitUtils::resolutionFunctionForVec->sigmaPt((recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol) )*std::pow(dmdpt[i],2), charge[i] );
280 
281  fullMassRes +=
282  ptVariance*std::pow(dmdpt[i],2) +
283  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
284  phiVariance*std::pow(dmdphi[i],2) +
285 
286  // These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
287  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
288  2*pt_phi*dmdpt[i]*dmdphi[i] +
289  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
290 
291  pt1_pt2*dmdpt[0]*dmdpt[1] +
292  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
293  phi1_phi2*dmdphi[0]*dmdphi[1] +
294 
295  // These are filled twice, because of the two combinations
296  pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
297  pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
298  pt12_phi21*dmdpt[0]*dmdphi[1] +
299  pt12_phi21*dmdpt[1]*dmdphi[0] +
300  cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
301  cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
302 
303  massRes1 += ptVariance*std::pow(dmdpt[i],2);
304  massRes2 += ptVariance*std::pow(dmdpt[i],2) +
305  cotgThetaVariance*std::pow(dmdcotgth[i],2);
306  massRes3 += ptVariance*std::pow(dmdpt[i],2) +
307  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
308  phiVariance*std::pow(dmdphi[i],2);
309  massRes4 += ptVariance*std::pow(dmdpt[i],2) +
310  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
311  phiVariance*std::pow(dmdphi[i],2) +
312  pt1_pt2*dmdpt[0]*dmdpt[1] +
313  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
314  2*pt_phi*dmdpt[i]*dmdphi[i] +
315  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i];
316  massRes5 += ptVariance*std::pow(dmdpt[i],2) +
317  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
318  phiVariance*std::pow(dmdphi[i],2) +
319  pt1_pt2*dmdpt[0]*dmdpt[1] +
320  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
321  2*pt_phi*dmdpt[i]*dmdphi[i] +
322  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
323  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
324  phi1_phi2*dmdphi[0]*dmdphi[1];
325  massRes6 += ptVariance*std::pow(dmdpt[i],2) +
326  cotgThetaVariance*std::pow(dmdcotgth[i],2) +
327  phiVariance*std::pow(dmdphi[i],2) +
328  pt1_pt2*dmdpt[0]*dmdpt[1] +
329  2*pt_cotgTheta*dmdpt[i]*dmdcotgth[i] +
330  2*pt_phi*dmdpt[i]*dmdphi[i] +
331  2*cotgTheta_phi*dmdcotgth[i]*dmdphi[i] +
332  cotgTheta1_cotgTheta2*dmdcotgth[0]*dmdcotgth[1] +
333  phi1_phi2*dmdphi[0]*dmdphi[1] +
334  pt12_cotgTheta21*dmdpt[0]*dmdcotgth[1] +
335  pt12_cotgTheta21*dmdpt[1]*dmdcotgth[0] +
336  pt12_phi21*dmdpt[0]*dmdphi[1] +
337  pt12_phi21*dmdpt[1]*dmdphi[0] +
338  cotgTheta12_phi21*dmdcotgth[0]*dmdphi[1] +
339  cotgTheta12_phi21*dmdcotgth[1]*dmdphi[0];
340 
341  massResPtAndPt12 += ptVariance*std::pow(dmdpt[i],2) + pt1_pt2*dmdpt[0]*dmdpt[1];
342 
343  // Derivatives
344  mapHisto_["DerivativePt"]->Fill( *(recMu[i]), dmdpt[i], charge[i] );
345  mapHisto_["DerivativeCotgTheta"]->Fill( *(recMu[i]), dmdcotgth[i], charge[i] );
346  mapHisto_["DerivativePhi"]->Fill( *(recMu[i]), dmdphi[i], charge[i] );
347  }
348  // Fill the complete resolution function with covariance terms
349  mapHisto_["FullMassResolution"]->Fill( *(recMu[0]), fullMassRes, charge[0]);
350  mapHisto_["FullMassResolution"]->Fill( *(recMu[1]), fullMassRes, charge[1]);
351 
352  mapHisto_["MassRes1"]->Fill( *(recMu[0]), massRes1, charge[0] );
353  mapHisto_["MassRes1"]->Fill( *(recMu[1]), massRes1, charge[1] );
354  mapHisto_["MassRes2"]->Fill( *(recMu[0]), massRes2, charge[0] );
355  mapHisto_["MassRes2"]->Fill( *(recMu[1]), massRes2, charge[1] );
356  mapHisto_["MassRes3"]->Fill( *(recMu[0]), massRes3, charge[0] );
357  mapHisto_["MassRes3"]->Fill( *(recMu[1]), massRes3, charge[1] );
358  mapHisto_["MassRes4"]->Fill( *(recMu[0]), massRes4, charge[0] );
359  mapHisto_["MassRes4"]->Fill( *(recMu[1]), massRes4, charge[1] );
360  mapHisto_["MassRes5"]->Fill( *(recMu[0]), massRes5, charge[0] );
361  mapHisto_["MassRes5"]->Fill( *(recMu[1]), massRes5, charge[1] );
362  mapHisto_["MassRes6"]->Fill( *(recMu[0]), massRes6, charge[0] );
363  mapHisto_["MassRes6"]->Fill( *(recMu[1]), massRes6, charge[1] );
364  mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[0]), massResPtAndPt12, charge[0] );
365  mapHisto_["MassResPtAndPt12"]->Fill( *(recMu[1]), massResPtAndPt12, charge[1] );
366  }
367  else {
368  // Fill the covariances histograms
369  mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
370  }
371  }
372  } // end if resonance
373  }
374 // else {
375 //
376 // // Loop on the recMuons
377 // std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
378 // for ( ; recMuon!=muons.end(); ++recMuon ) {
379 // int charge = recMuon->charge();
380 //
381 // lorentzVector recMu(recMuon->p4());
382 //
383 // // Find the matching MC muon
384 // const HepMC::GenEvent* Evt = evtMC->GetEvent();
385 // //Loop on generated particles
386 // std::map<double, lorentzVector> genAssocMap;
387 // HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
388 // for( ; part!=Evt->particles_end(); ++part ) {
389 // if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
390 // lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
391 // (*part)->momentum().pz(),(*part)->momentum().e()));
392 //
393 // double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
394 // ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
395 //
396 // // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
397 // // Only in this case we consider it matching.
398 // if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
399 // }
400 // }
401 // // Take the closest in deltaR
402 // lorentzVector genMu(genAssocMap.begin()->second);
403 //
404 // // Histograms with genParticles characteristics
405 // // --------------------------------------------
406 //
407 // if(checkDeltaR(genMu,recMu)){
408 // mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
409 // mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
410 // mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
411 // +cos(recMu.Theta())/sin(recMu.Theta())),charge);
412 // mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
413 // mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
414 // }
415 // }
416 }
417 
418 void ResolutionAnalyzer::fillHistoMap() {
419 
420  outputFile_->cd();
421 
422  // Resonances
423  // If no Z is required, use a smaller mass range.
424  double minMass = 0.;
425  double maxMass = 200.;
426  if( MuScleFitUtils::resfind[0] != 1 ) {
427  maxMass = 30.;
428  }
429  mapHisto_["GenMother"] = new HParticle(outputFile_, "GenMother", minMass, maxMass);
430  mapHisto_["SimResonance"] = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
431  mapHisto_["RecoResonance"] = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
432 
433  // Resonance muons
434  mapHisto_["GenMotherMuons"] = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
435  mapHisto_["SimResonanceMuons"] = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
436  mapHisto_["RecoResonanceMuons"] = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
437 
438  // Deltas between resonance muons
439  mapHisto_["DeltaGenMotherMuons"] = new HDelta (outputFile_, "DeltaGenMotherMuons");
440  mapHisto_["DeltaSimResonanceMuons"] = new HDelta (outputFile_, "DeltaSimResonanceMuons");
441  mapHisto_["DeltaRecoResonanceMuons"] = new HDelta (outputFile_, "DeltaRecoResonanceMuons");
442 
443  // //Reconstructed muon kinematics
444  // //-----------------------------
445  // mapHisto_["hRecBestMu"] = new HParticle ("hRecBestMu");
446  // mapHisto_["hRecBestMu_Acc"] = new HParticle ("hRecBestMu_Acc");
447  // mapHisto_["hDeltaRecBestMu"] = new HDelta ("hDeltaRecBestMu");
448 
449  // mapHisto_["hRecBestRes"] = new HParticle ("hRecBestRes");
450  // mapHisto_["hRecBestRes_Acc"] = new HParticle ("hRecBestRes_Acc");
451  // mapHisto_["hRecBestResVSMu"] = new HMassVSPart ("hRecBestResVSMu");
452 
453  //Resolution VS muon kinematic
454  //----------------------------
455  mapHisto_["PtResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionGenVSMu");
456  mapHisto_["PtResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PtResolutionSimVSMu");
457  mapHisto_["EtaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionGenVSMu");
458  mapHisto_["EtaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "EtaResolutionSimVSMu");
459  mapHisto_["ThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionGenVSMu");
460  mapHisto_["ThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "ThetaResolutionSimVSMu");
461  mapHisto_["CotgThetaResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
462  mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "CotgThetaResolutionSimVSMu");
463  mapHisto_["PhiResolutionGenVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
464  mapHisto_["PhiResolutionSimVSMu"] = new HResolutionVSPart (outputFile_, "PhiResolutionSimVSMu");
465 
466  // Covariances between muons kinematic quantities
467  // ----------------------------------------------
468  double ptMax = ptMax_;
469 
470  // Mass resolution
471  // ---------------
472  mapHisto_["MassResolution"] = new HMassResolutionVSPart (outputFile_,"MassResolution");
473 
474  // mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
475 
476  // Mass resolution vs (pt, eta) of the muons from MC
477  massResolutionVsPtEta_ = new HCovarianceVSxy ( "Mass", "Mass", 100, 0., ptMax, 60, -3, 3 );
478  // Mass resolution vs (pt, eta) of the muons from function
479  recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
480  recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
481  deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
482  deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
483 
484  // Muons resolutions from resolution functions
485  // -------------------------------------------
486  int totBinsY = 60;
487  mapHisto_["hFunctionResolMass"] = new HFunctionResolution (outputFile_, "hFunctionResolMass", ptMax, totBinsY);
488  mapHisto_["hFunctionResolPt"] = new HFunctionResolution (outputFile_, "hFunctionResolPt", ptMax, totBinsY);
489  mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution (outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
490  mapHisto_["hFunctionResolPhi"] = new HFunctionResolution (outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
491 
492  if( readCovariances_ ) {
493  // Covariances read from file. Used to compare the terms in the expression of mass error
494  mapHisto_["ReadCovariances"] = new HCovarianceVSParts ( theCovariancesRootFileName_, "Covariance" );
495 
496  // Variances
497  mapHisto_["MassResolutionPt"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPt", ptMax);
498  mapHisto_["MassResolutionCotgTheta"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassCotgTheta", ptMax);
499  mapHisto_["MassResolutionPhi"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPhi", ptMax);
500  // Covariances
501  mapHisto_["MassResolutionPt-CotgTheta"] = new HFunctionResolution(outputFile_,"functionResolMassPt-CotgTheta", ptMax, totBinsY);
502  mapHisto_["MassResolutionPt-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassPt-Phi", ptMax, totBinsY);
503  mapHisto_["MassResolutionCotgTheta-Phi"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta-Phi", ptMax, totBinsY);
504  mapHisto_["MassResolutionPt1-Pt2"] = new HFunctionResolution(outputFile_,"functionResolMassPt1-Pt2", ptMax, totBinsY);
505  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
506  mapHisto_["MassResolutionPhi1-Phi2"] = new HFunctionResolution(outputFile_,"functionResolMassPhi1-Phi2", ptMax, totBinsY);
507  mapHisto_["MassResolutionPt12-CotgTheta21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
508  mapHisto_["MassResolutionPt12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassPt12-Phi21", ptMax, totBinsY);
509  mapHisto_["MassResolutionCotgTheta12-Phi21"] = new HFunctionResolution(outputFile_,"functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
510 
511  mapHisto_["sigmaPtFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPtFromVariance", ptMax, totBinsY);
512  mapHisto_["sigmaCotgThetaFromVariance"] = new HFunctionResolution(outputFile_,"sigmaCotgThetaFromVariance", ptMax, totBinsY);
513  mapHisto_["sigmaPhiFromVariance"] = new HFunctionResolution(outputFile_,"sigmaPhiFromVariance", ptMax, totBinsY);
514 
515  // Derivatives
516  mapHisto_["DerivativePt"] = new HFunctionResolution(outputFile_,"derivativePt", ptMax);
517  mapHisto_["DerivativeCotgTheta"] = new HFunctionResolution(outputFile_,"derivativeCotgTheta", ptMax);
518  mapHisto_["DerivativePhi"] = new HFunctionResolution(outputFile_,"derivativePhi", ptMax);
519 
520  // Pt term from function
521  mapHisto_["MassResolutionPtFromFunction"] = new HFunctionResolutionVarianceCheck(outputFile_,"functionResolMassPtFromFunction", ptMax);
522 
523  mapHisto_["FullMassResolution"] = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
524  mapHisto_["MassRes1"] = new HFunctionResolution(outputFile_, "massRes1", ptMax);
525  mapHisto_["MassRes2"] = new HFunctionResolution(outputFile_, "massRes2", ptMax);
526  mapHisto_["MassRes3"] = new HFunctionResolution(outputFile_, "massRes3", ptMax);
527  mapHisto_["MassRes4"] = new HFunctionResolution(outputFile_, "massRes4", ptMax);
528  mapHisto_["MassRes5"] = new HFunctionResolution(outputFile_, "massRes5", ptMax);
529  mapHisto_["MassRes6"] = new HFunctionResolution(outputFile_, "massRes6", ptMax);
530  mapHisto_["MassResPtAndPt12"] = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
531  }
532  else {
533  mapHisto_["Covariances"] = new HCovarianceVSParts ( outputFile_, "Covariance", ptMax );
534  }
535 }
536 
537 void ResolutionAnalyzer::writeHistoMap() {
538  for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin();
539  histo!=mapHisto_.end(); histo++) {
540  (*histo).second->Write();
541  }
542  outputFile_->cd();
543  massResolutionVsPtEta_->Write();
544  recoPtVsgenPt_->Write();
545  recoPtVsgenPtEta12_->Write();
546  deltaPtOverPt_->Write();
547  deltaPtOverPtForEta12_->Write();
548 }
549 
550 bool ResolutionAnalyzer::checkDeltaR(const reco::Particle::LorentzVector & genMu, const reco::Particle::LorentzVector & recMu){
551  //first is always mu-, second is always mu+
552  double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) +
553  ((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
554  if(deltaR<0.01)
555  return true;
556  else if (debug_>0)
557  std::cout<<"Reco muon "<<recMu<<" with eta "<<recMu.Eta()<<" and phi "<<recMu.Phi()<<std::endl
558  <<" DOES NOT MATCH with generated muon from resonance: "<<std::endl
559  <<genMu<<" with eta "<<genMu.Eta()<<" and phi "<<genMu.Phi()<<std::endl;
560  return false;
561 }
562 
563 //define this as a plug-in
564 DEFINE_FWK_MODULE(ResolutionAnalyzer);
565 
566 #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:96
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
ResolutionAnalyzer::mapHisto_
std::map< std::string, Histograms * > mapHisto_
Definition: ResolutionAnalyzer.h:110
ResolutionAnalyzer::theRootFileName_
std::string theRootFileName_
Definition: ResolutionAnalyzer.h:107
jptDQMConfig_cff.ptMax
int ptMax
Definition: jptDQMConfig_cff.py:25
MuScleFitUtils::parResol
static std::vector< double > parResol
Definition: MuScleFitUtils.h:170
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
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
trackerHits.histo
tuple histo
Definition: trackerHits.py:12
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:8
HCovarianceVSxy::Fill
virtual void Fill(const double &x, const double &y, const double &a, const double &b)
Definition: Histograms.h:1611
ResolutionAnalyzer::massResolutionVsPtEta_
HCovarianceVSxy * massResolutionVsPtEta_
Definition: ResolutionAnalyzer.h:122
ResolutionAnalyzer::fillHistoMap
void fillHistoMap()
Used to fill the map with the histograms needed.
Definition: ResolutionAnalyzer.cc:418
iEvent
int iEvent
Definition: GenABIO.cc:243
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:116
MuScleFitUtils::massResolution
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
ResolutionAnalyzer::writeHistoMap
void writeHistoMap()
Writes the histograms in the map.
Definition: ResolutionAnalyzer.cc:537
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
MuonPairVector
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
Definition: RootTreeHandler.h:12
ResolutionAnalyzer::theCovariancesRootFileName_
std::string theCovariancesRootFileName_
Definition: ResolutionAnalyzer.h:108
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:550
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:130
MuScleFitUtils::mMu2
static const double mMu2
Definition: MuScleFitUtils.h:117
resolutionFunctionService
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:88
MuScleFitUtils::deltaPhi
static double deltaPhi(const double &phi1, const double &phi2)
Definition: MuScleFitUtils.h:88
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
MuScleFitUtils::resfind
static std::vector< int > resfind
Definition: MuScleFitUtils.h:182
MuScleFitUtils::resolutionFunctionForVec
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
Definition: MuScleFitUtils.h:132
reco::Particle::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:26
HResolutionVSPart
A set of histograms for resolution.
Definition: Histograms.h:857
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:131
ResolutionAnalyzer::ResolutionAnalyzer
ResolutionAnalyzer(const edm::ParameterSet &)
Definition: ResolutionAnalyzer.cc:19
Generated for CMSSW Reference Manual by doxygen 1.8.5