1 #ifndef RESOLUTIONANALYZER_CC
2 #define RESOLUTIONANALYZER_CC
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")) {
34 int resolFitType = iConfig.
getParameter<
int>(
"ResolFitType");
72 typedef std::vector<std::pair<lorentzVector, lorentzVector> >
MuonPairVector;
76 std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
78 MuonPairVector::iterator savedPair = savedPairVector.begin();
79 MuonPairVector::iterator genPair = genPairVector.begin();
80 std::cout <<
"Starting loop on " << savedPairVector.size() <<
" muons" << std::endl;
81 for (; savedPair != savedPairVector.end(); ++savedPair, ++genPair) {
98 mapHisto_[
"DeltaGenMotherMuons"]->Fill(genPair->first, genPair->second);
99 mapHisto_[
"GenMotherMuons"]->Fill(genPair->first);
100 mapHisto_[
"GenMotherMuons"]->Fill(genPair->second);
105 mapHisto_[
"PtResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Pt() + recMu1.Pt()) / genPair->first.Pt(), -1);
106 mapHisto_[
"ThetaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Theta() + recMu1.Theta()), -1);
107 mapHisto_[
"CotgThetaResolutionGenVSMu"]->Fill(
109 (-
cos(genPair->first.Theta()) /
sin(genPair->first.Theta()) +
cos(recMu1.Theta()) /
sin(recMu1.Theta())),
111 mapHisto_[
"EtaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Eta() + recMu1.Eta()), -1);
116 deltaPtOverPt_->Fill((recMu1.Pt() - genPair->first.Pt()) / genPair->first.Pt());
117 if (fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2) {
124 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(
128 (-
cos(genPair->second.Theta()) /
sin(genPair->second.Theta()) +
cos(recMu2.Theta()) /
sin(recMu2.Theta())),
130 mapHisto_[
"EtaResolutionGenVSMu"]->Fill(recMu2, (-genPair->second.Eta() + recMu2.Eta()), +1);
135 deltaPtOverPt_->Fill((recMu2.Pt() - genPair->second.Pt()) / genPair->second.Pt());
136 if (fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2) {
146 if (genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
148 double recoMass = (recMu1 + recMu2).
mass();
149 double genMass = (genPair->first + genPair->second).
mass();
151 mapHisto_[
"MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
155 mapHisto_[
"RecoResonance"]->Fill(recoResonance);
156 mapHisto_[
"DeltaRecoResonanceMuons"]->Fill(recMu1, recMu2);
157 mapHisto_[
"RecoResonanceMuons"]->Fill(recMu1);
158 mapHisto_[
"RecoResonanceMuons"]->Fill(recMu2);
163 double diffMass = recoMass - genMass;
165 double pt1 = recMu1.pt();
166 double eta1 = recMu1.eta();
167 double pt2 = recMu2.pt();
168 double eta2 = recMu2.eta();
170 if (diffMass == diffMass) {
174 std::cout <<
"Error, there is a nan: recoMass = " << recoMass <<
", genMass = " << genMass << std::endl;
189 mapHisto_[
"hFunctionResolCotgTheta"]->Fill(
201 mapHisto_[
"hFunctionResolCotgTheta"]->Fill(
213 double mass = (recMu1 + recMu2).
mass();
214 double pt1 = recMu1.Pt();
215 double phi1 = recMu1.Phi();
216 double eta1 = recMu1.Eta();
217 double theta1 = 2 * atan(
exp(-
eta1));
218 double pt2 = recMu2.Pt();
219 double phi2 = recMu2.Phi();
220 double eta2 = recMu2.Eta();
221 double theta2 = 2 * atan(
exp(-
eta2));
247 double dmdpt[2] = {dmdpt1 * recMu1.Pt(), dmdpt2 * recMu2.Pt()};
248 double dmdphi[2] = {dmdphi1, dmdphi2};
249 double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
256 double fullMassRes = 0.;
257 double massRes1 = 0.;
258 double massRes2 = 0.;
259 double massRes3 = 0.;
260 double massRes4 = 0.;
261 double massRes5 = 0.;
262 double massRes6 = 0.;
263 double massResPtAndPt12 = 0.;
265 for (
int i = 0;
i < 2; ++
i) {
266 double ptVariance =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt");
267 double cotgThetaVariance =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"CotgTheta");
268 double phiVariance =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Phi");
269 double pt_cotgTheta =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt-CotgTheta");
270 double pt_phi =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt-Phi");
271 double cotgTheta_phi =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"CotgTheta-Phi");
273 double pt1_pt2 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt1-Pt2");
274 double cotgTheta1_cotgTheta2 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"CotgTheta1-CotgTheta2");
275 double phi1_phi2 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Phi1-Phi2");
276 double pt12_cotgTheta21 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt12-CotgTheta21");
277 double pt12_phi21 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"Pt12-Phi21");
278 double cotgTheta12_phi21 =
mapHisto_[
"ReadCovariances"]->Get(*(recMu[
i]),
"CotgTheta12-Phi21");
282 mapHisto_[
"MassResolutionCotgTheta"]->Fill(
285 mapHisto_[
"MassResolutionPt-CotgTheta"]->Fill(
286 *(recMu[
i]), 2 * pt_cotgTheta * dmdpt[
i] * dmdcotgth[
i],
charge[
i]);
287 mapHisto_[
"MassResolutionPt-Phi"]->Fill(*(recMu[
i]), 2 * pt_phi * dmdpt[
i] * dmdphi[
i],
charge[
i]);
288 mapHisto_[
"MassResolutionCotgTheta-Phi"]->Fill(
289 *(recMu[
i]), 2 * cotgTheta_phi * dmdcotgth[
i] * dmdphi[
i],
charge[
i]);
291 mapHisto_[
"MassResolutionPt1-Pt2"]->Fill(*(recMu[
i]), pt1_pt2 * dmdpt[0] * dmdpt[1],
charge[
i]);
292 mapHisto_[
"MassResolutionCotgTheta1-CotgTheta2"]->Fill(
293 *(recMu[
i]), cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1],
charge[
i]);
294 mapHisto_[
"MassResolutionPhi1-Phi2"]->Fill(*(recMu[
i]), phi1_phi2 * dmdphi[0] * dmdphi[1],
charge[
i]);
296 mapHisto_[
"MassResolutionPt12-CotgTheta21"]->Fill(
297 *(recMu[
i]), pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1],
charge[
i]);
298 mapHisto_[
"MassResolutionPt12-CotgTheta21"]->Fill(
299 *(recMu[
i]), pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0],
charge[
i]);
300 mapHisto_[
"MassResolutionPt12-Phi21"]->Fill(*(recMu[
i]), pt12_phi21 * dmdpt[0] * dmdphi[1],
charge[
i]);
301 mapHisto_[
"MassResolutionPt12-Phi21"]->Fill(*(recMu[
i]), pt12_phi21 * dmdpt[1] * dmdphi[0],
charge[
i]);
302 mapHisto_[
"MassResolutionCotgTheta12-Phi21"]->Fill(
303 *(recMu[
i]), cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1],
charge[
i]);
304 mapHisto_[
"MassResolutionCotgTheta12-Phi21"]->Fill(
305 *(recMu[
i]), cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0],
charge[
i]);
313 mapHisto_[
"MassResolutionPtFromFunction"]->Fill(
320 fullMassRes += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2) +
324 2 * pt_cotgTheta * dmdpt[
i] * dmdcotgth[
i] + 2 * pt_phi * dmdpt[
i] * dmdphi[
i] +
325 2 * cotgTheta_phi * dmdcotgth[
i] * dmdphi[
i] +
327 pt1_pt2 * dmdpt[0] * dmdpt[1] + cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] +
328 phi1_phi2 * dmdphi[0] * dmdphi[1] +
331 pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
332 pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
333 cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
335 massRes1 += ptVariance *
std::pow(dmdpt[
i], 2);
336 massRes2 += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2);
337 massRes3 += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2) +
339 massRes4 += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2) +
340 phiVariance *
std::pow(dmdphi[
i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
341 2 * pt_cotgTheta * dmdpt[
i] * dmdcotgth[
i] + 2 * pt_phi * dmdpt[
i] * dmdphi[
i] +
342 2 * cotgTheta_phi * dmdcotgth[
i] * dmdphi[
i];
343 massRes5 += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2) +
344 phiVariance *
std::pow(dmdphi[
i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
345 2 * pt_cotgTheta * dmdpt[
i] * dmdcotgth[
i] + 2 * pt_phi * dmdpt[
i] * dmdphi[
i] +
346 2 * cotgTheta_phi * dmdcotgth[
i] * dmdphi[
i] +
347 cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1];
348 massRes6 += ptVariance *
std::pow(dmdpt[
i], 2) + cotgThetaVariance *
std::pow(dmdcotgth[
i], 2) +
349 phiVariance *
std::pow(dmdphi[
i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
350 2 * pt_cotgTheta * dmdpt[
i] * dmdcotgth[
i] + 2 * pt_phi * dmdpt[
i] * dmdphi[
i] +
351 2 * cotgTheta_phi * dmdcotgth[
i] * dmdphi[
i] +
352 cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1] +
353 pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
354 pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
355 cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
357 massResPtAndPt12 += ptVariance *
std::pow(dmdpt[
i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1];
365 mapHisto_[
"FullMassResolution"]->Fill(*(recMu[0]), fullMassRes,
charge[0]);
366 mapHisto_[
"FullMassResolution"]->Fill(*(recMu[1]), fullMassRes,
charge[1]);
380 mapHisto_[
"MassResPtAndPt12"]->Fill(*(recMu[0]), massResPtAndPt12,
charge[0]);
381 mapHisto_[
"MassResPtAndPt12"]->Fill(*(recMu[1]), massResPtAndPt12,
charge[1]);
384 mapHisto_[
"Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
475 mapHisto_[
"CotgThetaResolutionGenVSMu"] =
497 deltaPtOverPt_ =
new TH1D(
"deltaPtOverPt",
"deltaPtOverPt", 100, -0.1, 0.1);
519 mapHisto_[
"MassResolutionPt-CotgTheta"] =
523 mapHisto_[
"MassResolutionCotgTheta-Phi"] =
527 mapHisto_[
"MassResolutionCotgTheta1-CotgTheta2"] =
531 mapHisto_[
"MassResolutionPt12-CotgTheta21"] =
535 mapHisto_[
"MassResolutionCotgTheta12-Phi21"] =
539 mapHisto_[
"sigmaCotgThetaFromVariance"] =
549 mapHisto_[
"MassResolutionPtFromFunction"] =
568 (*histo).second->Write();
583 ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
587 std::cout <<
"Reco muon " << recMu <<
" with eta " << recMu.Eta() <<
" and phi " << recMu.Phi() << std::endl
588 <<
" DOES NOT MATCH with generated muon from resonance: " << std::endl
589 << genMu <<
" with eta " << genMu.Eta() <<
" and phi " << genMu.Phi() << std::endl;
596 #endif // RESOLUTIONANALYZER_CC