File indexing completed on 2021-02-14 13:31:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107 #include "FWCore/Framework/interface/EDLooper.h"
0108 #include "FWCore/Framework/interface/LooperFactory.h"
0109 #include "FWCore/Utilities/interface/InputTag.h"
0110 #include "FWCore/Framework/interface/Frameworkfwd.h"
0111 #include "FWCore/Framework/interface/Event.h"
0112 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0113
0114 #include <CLHEP/Vector/LorentzVector.h>
0115 #include <memory>
0116
0117 #include <vector>
0118
0119 #include "FWCore/Framework/interface/EventSetup.h"
0120 #include "FWCore/Framework/interface/ESHandle.h"
0121 #include "FWCore/ServiceRegistry/interface/Service.h"
0122
0123 #include "MuScleFitBase.h"
0124 #include "Histograms.h"
0125 #include "MuScleFitPlotter.h"
0126 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
0127 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
0128 #include "MuonAnalysis/MomentumScaleCalibration/interface/Muon.h"
0129 #include "MuonAnalysis/MomentumScaleCalibration/interface/Event.h"
0130 #include "MuScleFitMuonSelector.h"
0131
0132 #include "DataFormats/TrackReco/interface/Track.h"
0133 #include "DataFormats/MuonReco/interface/Muon.h"
0134 #include "DataFormats/MuonReco/interface/CaloMuon.h"
0135 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0136
0137 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0138 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0139
0140 #include "DataFormats/PatCandidates/interface/Muon.h"
0141 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0142 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0143 #include "DataFormats/Common/interface/TriggerResults.h"
0144 #include "FWCore/Common/interface/TriggerNames.h"
0145
0146 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0147
0148 #include "HepPDT/defs.h"
0149 #include "HepPDT/TableBuilder.hh"
0150 #include "HepPDT/ParticleDataTable.hh"
0151
0152 #include "HepMC/GenParticle.h"
0153 #include "HepMC/GenEvent.h"
0154
0155 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0156 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0157 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0158 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0159
0160 #include "DataFormats/VertexReco/interface/Vertex.h"
0161 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0162 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0163
0164 #include "TFile.h"
0165 #include "TTree.h"
0166 #include "TMinuit.h"
0167
0168
0169
0170
0171 #ifdef USE_CALLGRIND
0172 #include "valgrind/callgrind.h"
0173 #endif
0174
0175
0176
0177
0178
0179 namespace edm {
0180 class ParameterSet;
0181 class Event;
0182 class EventSetup;
0183 }
0184
0185 class MuScleFit : public edm::EDLooper, MuScleFitBase {
0186 public:
0187
0188
0189 MuScleFit(const edm::ParameterSet& pset);
0190
0191
0192
0193 ~MuScleFit() override;
0194
0195
0196
0197 void beginOfJobInConstructor();
0198
0199
0200 void endOfJob() override;
0201
0202 void startingNewLoop(unsigned int iLoop) override;
0203
0204 edm::EDLooper::Status endOfLoop(const edm::EventSetup& eventSetup, unsigned int iLoop) override;
0205 virtual void endOfFastLoop(const unsigned int iLoop);
0206
0207 edm::EDLooper::Status duringLoop(const edm::Event& event, const edm::EventSetup& eventSetup) override;
0208
0209
0210
0211
0212 virtual void duringFastLoop();
0213
0214 template <typename T>
0215 std::vector<MuScleFitMuon> fillMuonCollection(const std::vector<T>& tracks);
0216
0217 private:
0218 protected:
0219
0220
0221
0222
0223 void selectMuons(const edm::Event& event);
0224
0225
0226
0227
0228
0229 void selectMuons(const int maxEvents, const TString& treeFileName);
0230
0231
0232 template <typename T>
0233 void takeSelectedMuonType(const T& muon, std::vector<reco::Track>& tracks);
0234
0235 bool selGlobalMuon(const pat::Muon* aMuon);
0236 bool selTrackerMuon(const pat::Muon* aMuon);
0237
0238
0239 bool checkDeltaR(reco::Particle::LorentzVector& genMu, reco::Particle::LorentzVector& recMu);
0240
0241 void fillComparisonHistograms(const reco::Particle::LorentzVector& genMu,
0242 const reco::Particle::LorentzVector& recoMu,
0243 const std::string& inputName,
0244 const int charge);
0245
0246
0247 void applySmearing(reco::Particle::LorentzVector& mu);
0248
0249 void applyBias(reco::Particle::LorentzVector& mu, const int charge);
0250
0251
0252
0253
0254
0255 void checkParameters();
0256
0257 MuonServiceProxy* theService;
0258
0259
0260
0261 int numberOfSimTracks;
0262 int numberOfSimMuons;
0263 int numberOfSimVertices;
0264 int numberOfEwkZ;
0265
0266 bool ifHepMC;
0267 bool ifGenPart;
0268
0269
0270
0271 double minResMass_hwindow[6];
0272 double maxResMass_hwindow[6];
0273
0274
0275
0276 unsigned int maxLoopNumber;
0277 unsigned int loopCounter;
0278
0279 bool fastLoop;
0280
0281 MuScleFitPlotter* plotter;
0282
0283
0284
0285 reco::Particle::LorentzVector recMu1, recMu2;
0286 MuScleFitMuon recMuScleMu1, recMuScleMu2;
0287 int iev;
0288 int totalEvents_;
0289
0290 bool compareToSimTracks_;
0291 edm::InputTag simTracksCollection_;
0292 bool PATmuons_;
0293 std::string genParticlesName_;
0294
0295
0296 std::string inputRootTreeFileName_;
0297
0298 std::string outputRootTreeFileName_;
0299
0300 int maxEventsFromRootTree_;
0301
0302 std::string triggerResultsLabel_;
0303 std::string triggerResultsProcess_;
0304 std::vector<std::string> triggerPath_;
0305 bool negateTrigger_;
0306 bool saveAllToTree_;
0307
0308
0309 edm::InputTag puInfoSrc_;
0310 edm::InputTag vertexSrc_;
0311
0312 std::unique_ptr<MuScleFitMuonSelector> muonSelector_;
0313 };
0314
0315 template <typename T>
0316 std::vector<MuScleFitMuon> MuScleFit::fillMuonCollection(const std::vector<T>& tracks) {
0317 std::vector<MuScleFitMuon> muons;
0318 typename std::vector<T>::const_iterator track;
0319 for (track = tracks.begin(); track != tracks.end(); ++track) {
0320 reco::Particle::LorentzVector mu;
0321 mu = reco::Particle::LorentzVector(
0322 track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
0323
0324
0325 MuScleFitUtils::goodmuon++;
0326 if (debug_ > 0)
0327 std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon << ": initial value Pt = " << mu.Pt()
0328 << std::endl;
0329
0330 applySmearing(mu);
0331 applyBias(mu, track->charge());
0332 if (debug_ > 0)
0333 std::cout << "track charge: " << track->charge() << std::endl;
0334
0335 Double_t hitsTk = track->innerTrack()->hitPattern().numberOfValidTrackerHits();
0336 Double_t hitsMuon = track->innerTrack()->hitPattern().numberOfValidMuonHits();
0337 Double_t ptError = track->innerTrack()->ptError();
0338 MuScleFitMuon muon(mu, track->charge(), ptError, hitsTk, hitsMuon, false);
0339 if (debug_ > 0) {
0340 std::cout << "[MuScleFit::fillMuonCollection]" << std::endl;
0341 std::cout << " muon = " << muon << std::endl;
0342 }
0343
0344
0345
0346 muons.push_back(muon);
0347 }
0348 return muons;
0349 }
0350
0351 template <typename T>
0352 void MuScleFit::takeSelectedMuonType(const T& muon, std::vector<reco::Track>& tracks) {
0353
0354
0355
0356 if (muon->isGlobalMuon() && theMuonType_ == 1)
0357 tracks.push_back(*(muon->globalTrack()));
0358 else if (muon->isStandAloneMuon() && theMuonType_ == 2)
0359 tracks.push_back(*(muon->outerTrack()));
0360 else if (muon->isTrackerMuon() && theMuonType_ == 3)
0361 tracks.push_back(*(muon->innerTrack()));
0362
0363 else if (theMuonType_ == 10 && !(muon->isStandAloneMuon()))
0364 tracks.push_back(*(muon->innerTrack()));
0365 else if (theMuonType_ == 11 && muon->isGlobalMuon())
0366 tracks.push_back(*(muon->innerTrack()));
0367 else if (theMuonType_ == 13 && muon->isTrackerMuon())
0368 tracks.push_back(*(muon->innerTrack()));
0369 }
0370
0371
0372
0373 MuScleFit::MuScleFit(const edm::ParameterSet& pset) : MuScleFitBase(pset), totalEvents_(0) {
0374 MuScleFitUtils::debug = debug_;
0375 if (debug_ > 0)
0376 std::cout << "[MuScleFit]: Constructor" << std::endl;
0377
0378 if ((theMuonType_ < -4 || theMuonType_ > 5) && theMuonType_ < 10) {
0379 std::cout << "[MuScleFit]: Unknown muon type! Aborting." << std::endl;
0380 abort();
0381 }
0382
0383 loopCounter = 0;
0384
0385
0386
0387 minResMass_hwindow[0] = 71.1876;
0388 maxResMass_hwindow[0] = 111.188;
0389 minResMass_hwindow[1] = 10.15;
0390 maxResMass_hwindow[1] = 10.55;
0391 minResMass_hwindow[2] = 9.8;
0392 maxResMass_hwindow[2] = 10.2;
0393 minResMass_hwindow[3] = 9.25;
0394 maxResMass_hwindow[3] = 9.65;
0395 minResMass_hwindow[4] = 3.58;
0396 maxResMass_hwindow[4] = 3.78;
0397 minResMass_hwindow[5] = 3.0;
0398 maxResMass_hwindow[5] = 3.2;
0399
0400
0401
0402 maxLoopNumber = pset.getUntrackedParameter<int>("maxLoopNumber", 2);
0403 fastLoop = pset.getUntrackedParameter<bool>("FastLoop", true);
0404
0405
0406 MuScleFitUtils::doResolFit = pset.getParameter<std::vector<int> >("doResolFit");
0407 MuScleFitUtils::doScaleFit = pset.getParameter<std::vector<int> >("doScaleFit");
0408 MuScleFitUtils::doCrossSectionFit = pset.getParameter<std::vector<int> >("doCrossSectionFit");
0409 MuScleFitUtils::doBackgroundFit = pset.getParameter<std::vector<int> >("doBackgroundFit");
0410
0411
0412
0413 int biasType = pset.getParameter<int>("BiasType");
0414 MuScleFitUtils::BiasType = biasType;
0415
0416 MuScleFitUtils::biasFunction = scaleFunctionVecService(biasType);
0417 int smearType = pset.getParameter<int>("SmearType");
0418 MuScleFitUtils::SmearType = smearType;
0419 MuScleFitUtils::smearFunction = smearFunctionService(smearType);
0420
0421
0422
0423 int resolFitType = pset.getParameter<int>("ResolFitType");
0424 MuScleFitUtils::ResolFitType = resolFitType;
0425 MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType);
0426 MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolFitType);
0427 int scaleType = pset.getParameter<int>("ScaleFitType");
0428 MuScleFitUtils::ScaleFitType = scaleType;
0429 MuScleFitUtils::scaleFunction = scaleFunctionService(scaleType);
0430 MuScleFitUtils::scaleFunctionForVec = scaleFunctionVecService(scaleType);
0431
0432
0433
0434 MuScleFitUtils::parBias = pset.getParameter<std::vector<double> >("parBias");
0435 MuScleFitUtils::parSmear = pset.getParameter<std::vector<double> >("parSmear");
0436 MuScleFitUtils::parResol = pset.getParameter<std::vector<double> >("parResol");
0437 MuScleFitUtils::parResolStep =
0438 pset.getUntrackedParameter<std::vector<double> >("parResolStep", std::vector<double>());
0439 MuScleFitUtils::parResolMin = pset.getUntrackedParameter<std::vector<double> >("parResolMin", std::vector<double>());
0440 MuScleFitUtils::parResolMax = pset.getUntrackedParameter<std::vector<double> >("parResolMax", std::vector<double>());
0441 MuScleFitUtils::parScale = pset.getParameter<std::vector<double> >("parScale");
0442 MuScleFitUtils::parScaleStep =
0443 pset.getUntrackedParameter<std::vector<double> >("parScaleStep", std::vector<double>());
0444 MuScleFitUtils::parScaleMin = pset.getUntrackedParameter<std::vector<double> >("parScaleMin", std::vector<double>());
0445 MuScleFitUtils::parScaleMax = pset.getUntrackedParameter<std::vector<double> >("parScaleMax", std::vector<double>());
0446 MuScleFitUtils::parCrossSection = pset.getParameter<std::vector<double> >("parCrossSection");
0447 MuScleFitUtils::parBgr = pset.getParameter<std::vector<double> >("parBgr");
0448 MuScleFitUtils::parResolFix = pset.getParameter<std::vector<int> >("parResolFix");
0449 MuScleFitUtils::parScaleFix = pset.getParameter<std::vector<int> >("parScaleFix");
0450 MuScleFitUtils::parCrossSectionFix = pset.getParameter<std::vector<int> >("parCrossSectionFix");
0451 MuScleFitUtils::parBgrFix = pset.getParameter<std::vector<int> >("parBgrFix");
0452 MuScleFitUtils::parResolOrder = pset.getParameter<std::vector<int> >("parResolOrder");
0453 MuScleFitUtils::parScaleOrder = pset.getParameter<std::vector<int> >("parScaleOrder");
0454 MuScleFitUtils::parCrossSectionOrder = pset.getParameter<std::vector<int> >("parCrossSectionOrder");
0455 MuScleFitUtils::parBgrOrder = pset.getParameter<std::vector<int> >("parBgrOrder");
0456
0457 MuScleFitUtils::resfind = pset.getParameter<std::vector<int> >("resfind");
0458 MuScleFitUtils::FitStrategy = pset.getParameter<int>("FitStrategy");
0459
0460
0461
0462 MuScleFitUtils::speedup = pset.getParameter<bool>("speedup");
0463
0464
0465 compareToSimTracks_ = pset.getParameter<bool>("compareToSimTracks");
0466 simTracksCollection_ = pset.getUntrackedParameter<edm::InputTag>("SimTracksCollection", edm::InputTag("g4SimHits"));
0467
0468 triggerResultsLabel_ = pset.getUntrackedParameter<std::string>("TriggerResultsLabel");
0469 triggerResultsProcess_ = pset.getUntrackedParameter<std::string>("TriggerResultsProcess");
0470 triggerPath_ = pset.getUntrackedParameter<std::vector<std::string> >("TriggerPath");
0471 negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
0472 saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
0473
0474
0475 puInfoSrc_ = pset.getUntrackedParameter<edm::InputTag>("PileUpSummaryInfo");
0476 vertexSrc_ = pset.getUntrackedParameter<edm::InputTag>("PrimaryVertexCollection");
0477
0478 PATmuons_ = pset.getUntrackedParameter<bool>("PATmuons", false);
0479 genParticlesName_ = pset.getUntrackedParameter<std::string>("GenParticlesName", "genParticles");
0480
0481
0482
0483 MuScleFitUtils::useProbsFile_ = pset.getUntrackedParameter<bool>("UseProbsFile", true);
0484
0485
0486 MuScleFitUtils::sherpa_ = pset.getUntrackedParameter<bool>("Sherpa", false);
0487
0488 MuScleFitUtils::rapidityBinsForZ_ = pset.getUntrackedParameter<bool>("RapidityBinsForZ", true);
0489
0490
0491 MuScleFitUtils::separateRanges_ = pset.getUntrackedParameter<bool>("SeparateRanges", true);
0492 MuScleFitUtils::maxMuonPt_ = pset.getUntrackedParameter<double>("MaxMuonPt", 100000000.);
0493 MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
0494 MuScleFitUtils::minMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MinMuonEtaFirstRange", -6.);
0495 MuScleFitUtils::maxMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaFirstRange", 6.);
0496 MuScleFitUtils::minMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MinMuonEtaSecondRange", -100.);
0497 MuScleFitUtils::maxMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaSecondRange", 100.);
0498 MuScleFitUtils::deltaPhiMinCut_ = pset.getUntrackedParameter<double>("DeltaPhiMinCut", -100.);
0499 MuScleFitUtils::deltaPhiMaxCut_ = pset.getUntrackedParameter<double>("DeltaPhiMaxCut", 100.);
0500
0501 MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("DebugMassResol", false);
0502
0503
0504
0505
0506 checkParameters();
0507
0508
0509
0510 if (MuScleFitUtils::SmearType > 0) {
0511 std::cout << "[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
0512 TF1 G("G", "[0]*exp(-0.5*pow(x,2))", -5., 5.);
0513 double norm = 1 / sqrt(2 * TMath::Pi());
0514 G.SetParameter(0, norm);
0515 for (int i = 0; i < 10000; i++) {
0516 for (int j = 0; j < 7; j++) {
0517 MuScleFitUtils::x[j][i] = G.GetRandom();
0518 }
0519 }
0520 }
0521 MuScleFitUtils::goodmuon = 0;
0522
0523 if (theMuonType_ > 0 && theMuonType_ < 4) {
0524 MuScleFitUtils::MuonTypeForCheckMassWindow = theMuonType_ - 1;
0525 MuScleFitUtils::MuonType = theMuonType_ - 1;
0526 } else if (theMuonType_ == 0 || theMuonType_ == 4 || theMuonType_ == 5 || theMuonType_ >= 10 || theMuonType_ == -1 ||
0527 theMuonType_ == -2 || theMuonType_ == -3 || theMuonType_ == -4) {
0528 MuScleFitUtils::MuonTypeForCheckMassWindow = 2;
0529 MuScleFitUtils::MuonType = 2;
0530 } else {
0531 std::cout << "Wrong muon type " << theMuonType_ << std::endl;
0532 exit(1);
0533 }
0534
0535
0536 if (theMuonType_ == 2) {
0537 MuScleFitUtils::rapidityBinsForZ_ = false;
0538 }
0539
0540
0541
0542 MuScleFitUtils::massWindowHalfWidth[0][0] = 20.;
0543 MuScleFitUtils::massWindowHalfWidth[0][1] = 0.35;
0544 MuScleFitUtils::massWindowHalfWidth[0][2] = 0.35;
0545 MuScleFitUtils::massWindowHalfWidth[0][3] = 0.35;
0546 MuScleFitUtils::massWindowHalfWidth[0][4] = 0.2;
0547 MuScleFitUtils::massWindowHalfWidth[0][5] = 0.2;
0548 MuScleFitUtils::massWindowHalfWidth[1][0] = 20.;
0549 MuScleFitUtils::massWindowHalfWidth[1][1] = 0.35;
0550 MuScleFitUtils::massWindowHalfWidth[1][2] = 0.35;
0551 MuScleFitUtils::massWindowHalfWidth[1][3] = 0.35;
0552 MuScleFitUtils::massWindowHalfWidth[1][4] = 0.2;
0553 MuScleFitUtils::massWindowHalfWidth[1][5] = 0.2;
0554 MuScleFitUtils::massWindowHalfWidth[2][0] = 20.;
0555 MuScleFitUtils::massWindowHalfWidth[2][1] = 0.35;
0556 MuScleFitUtils::massWindowHalfWidth[2][2] = 0.35;
0557 MuScleFitUtils::massWindowHalfWidth[2][3] = 0.35;
0558 MuScleFitUtils::massWindowHalfWidth[2][4] = 0.2;
0559 MuScleFitUtils::massWindowHalfWidth[2][5] = 0.2;
0560
0561 muonSelector_ = std::make_unique<MuScleFitMuonSelector>(theMuonLabel_,
0562 theMuonType_,
0563 PATmuons_,
0564 MuScleFitUtils::resfind,
0565 MuScleFitUtils::speedup,
0566 genParticlesName_,
0567 compareToSimTracks_,
0568 simTracksCollection_,
0569 MuScleFitUtils::sherpa_,
0570 debug_);
0571
0572 MuScleFitUtils::backgroundHandler =
0573 new BackgroundHandler(pset.getParameter<std::vector<int> >("BgrFitType"),
0574 pset.getParameter<std::vector<double> >("LeftWindowBorder"),
0575 pset.getParameter<std::vector<double> >("RightWindowBorder"),
0576 MuScleFitUtils::ResMass,
0577 MuScleFitUtils::massWindowHalfWidth[MuScleFitUtils::MuonTypeForCheckMassWindow]);
0578
0579 MuScleFitUtils::crossSectionHandler =
0580 new CrossSectionHandler(MuScleFitUtils::parCrossSection, MuScleFitUtils::resfind);
0581
0582
0583
0584
0585 MuScleFitUtils::normalizeLikelihoodByEventNumber_ =
0586 pset.getUntrackedParameter<bool>("NormalizeLikelihoodByEventNumber", true);
0587 if (debug_ > 0)
0588 std::cout << "End of MuScleFit constructor" << std::endl;
0589
0590 inputRootTreeFileName_ = pset.getParameter<std::string>("InputRootTreeFileName");
0591 outputRootTreeFileName_ = pset.getParameter<std::string>("OutputRootTreeFileName");
0592 maxEventsFromRootTree_ = pset.getParameter<int>("MaxEventsFromRootTree");
0593
0594 MuScleFitUtils::startWithSimplex_ = pset.getParameter<bool>("StartWithSimplex");
0595 MuScleFitUtils::computeMinosErrors_ = pset.getParameter<bool>("ComputeMinosErrors");
0596 MuScleFitUtils::minimumShapePlots_ = pset.getParameter<bool>("MinimumShapePlots");
0597
0598 beginOfJobInConstructor();
0599 }
0600
0601
0602
0603 MuScleFit::~MuScleFit() {
0604 if (debug_ > 0)
0605 std::cout << "[MuScleFit]: Destructor" << std::endl;
0606 std::cout << "Total number of analyzed events = " << totalEvents_ << std::endl;
0607
0608 if (!(outputRootTreeFileName_.empty())) {
0609
0610 if (!(inputRootTreeFileName_.empty() && (int(MuScleFitUtils::SavedPair.size()) != totalEvents_))) {
0611 std::cout << "Saving muon pairs to root tree" << std::endl;
0612 RootTreeHandler rootTreeHandler;
0613 if (MuScleFitUtils::speedup) {
0614
0615 if (debug_ > 0) {
0616 std::vector<MuonPair>::const_iterator it = muonPairs_.begin();
0617 std::cout << "[MuScleFit::~MuScleFit] (Destructor)" << std::endl;
0618 for (; it < muonPairs_.end(); ++it) {
0619 std::cout << " Debugging pairs that are going to be written to file" << std::endl;
0620 std::cout << " muon1 = " << it->mu1 << std::endl;
0621 std::cout << " muon2 = " << it->mu2 << std::endl;
0622 }
0623 }
0624 rootTreeHandler.writeTree(outputRootTreeFileName_, &(muonPairs_), theMuonType_, nullptr, saveAllToTree_);
0625 } else {
0626
0627 rootTreeHandler.writeTree(
0628 outputRootTreeFileName_, &(muonPairs_), theMuonType_, &(genMuonPairs_), saveAllToTree_);
0629 }
0630 } else {
0631 std::cout << "ERROR: events in the vector = " << MuScleFitUtils::SavedPair.size()
0632 << " != totalEvents = " << totalEvents_ << std::endl;
0633 }
0634 }
0635 }
0636
0637
0638
0639 void MuScleFit::beginOfJobInConstructor()
0640
0641
0642 {
0643 if (debug_ > 0)
0644 std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0645
0646 if (MuScleFitUtils::useProbsFile_) {
0647 readProbabilityDistributionsFromFile();
0648 }
0649
0650 if (debug_ > 0)
0651 std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0652
0653
0654
0655 for (unsigned int i = 0; i < (maxLoopNumber); i++) {
0656 std::stringstream ss;
0657 ss << i;
0658 std::string rootFileName = ss.str() + "_" + theRootFileName_;
0659 if (theCompressionSettings_ > -1) {
0660 theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE", "", theCompressionSettings_));
0661 } else {
0662 theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE"));
0663 }
0664 }
0665 if (debug_ > 0)
0666 std::cout << "[MuScleFit]: Root file created" << std::endl;
0667
0668 std::cout << "creating plotter" << std::endl;
0669 plotter = new MuScleFitPlotter(theGenInfoRootFileName_);
0670 plotter->debug = debug_;
0671 }
0672
0673
0674
0675 void MuScleFit::endOfJob() {
0676 if (debug_ > 0)
0677 std::cout << "[MuScleFit]: endOfJob" << std::endl;
0678 }
0679
0680
0681
0682 void MuScleFit::startingNewLoop(unsigned int iLoop) {
0683 if (debug_ > 0)
0684 std::cout << "[MuScleFit]: Starting loop # " << iLoop << std::endl;
0685
0686
0687
0688 MuScleFitUtils::goodmuon = 0;
0689
0690
0691
0692 MuScleFitUtils::counter_resprob = 0;
0693
0694
0695
0696 fillHistoMap(theFiles_[iLoop], iLoop);
0697
0698 loopCounter = iLoop;
0699 MuScleFitUtils::loopCounter = loopCounter;
0700
0701 iev = 0;
0702 MuScleFitUtils::iev_ = 0;
0703
0704 MuScleFitUtils::oldNormalization_ = 0;
0705 }
0706
0707
0708
0709 edm::EDLooper::Status MuScleFit::endOfLoop(const edm::EventSetup& eventSetup, unsigned int iLoop) {
0710 unsigned int iFastLoop = 1;
0711
0712
0713 if (!(inputRootTreeFileName_.empty())) {
0714 selectMuons(maxEventsFromRootTree_, inputRootTreeFileName_);
0715
0716 totalEvents_ = MuScleFitUtils::SavedPair.size();
0717 iFastLoop = 0;
0718 } else {
0719 endOfFastLoop(iLoop);
0720 }
0721
0722
0723 if (fastLoop == true) {
0724 for (; iFastLoop < maxLoopNumber; ++iFastLoop) {
0725 std::cout << "Starting fast loop number " << iFastLoop << std::endl;
0726
0727
0728
0729 startingNewLoop(iFastLoop);
0730
0731
0732
0733
0734 while (iev < totalEvents_) {
0735 if (iev % 50000 == 0) {
0736 std::cout << "Fast looping on event number " << iev << std::endl;
0737 }
0738
0739 duringFastLoop();
0740 }
0741 std::cout << "End of fast loop number " << iFastLoop << ". Ran on " << iev << " events" << std::endl;
0742 endOfFastLoop(iFastLoop);
0743 }
0744 }
0745
0746 if (iFastLoop >= maxLoopNumber - 1) {
0747 return kStop;
0748 } else {
0749 return kContinue;
0750 }
0751 }
0752
0753 void MuScleFit::endOfFastLoop(const unsigned int iLoop) {
0754
0755
0756 if (loopCounter == 0) {
0757
0758
0759 delete plotter;
0760 }
0761
0762 std::cout << "Ending loop # " << iLoop << std::endl;
0763
0764
0765
0766
0767 writeHistoMap(iLoop);
0768
0769
0770
0771
0772 TDirectory* likelihoodDir = theFiles_[iLoop]->mkdir("likelihood");
0773 likelihoodDir->cd();
0774 MuScleFitUtils::minimizeLikelihood();
0775
0776
0777 theFiles_[iLoop]->Close();
0778
0779 delete theFiles_[iLoop];
0780
0781
0782
0783 clearHistoMap();
0784 }
0785
0786
0787
0788 edm::EDLooper::Status MuScleFit::duringLoop(const edm::Event& event, const edm::EventSetup& eventSetup) {
0789 edm::Handle<edm::TriggerResults> triggerResults;
0790 event.getByLabel(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()), triggerResults);
0791
0792 bool isFired = false;
0793
0794 if (triggerPath_[0].empty())
0795 isFired = true;
0796 else if (triggerPath_[0] == "All") {
0797 isFired = triggerResults->accept();
0798 if (debug_ > 0)
0799 std::cout << "Trigger " << isFired << std::endl;
0800 } else {
0801 bool changed;
0802 HLTConfigProvider hltConfig;
0803 hltConfig.init(event.getRun(), eventSetup, triggerResultsProcess_, changed);
0804
0805 const edm::TriggerNames& triggerNames = event.triggerNames(*triggerResults);
0806
0807 for (unsigned i = 0; i < triggerNames.size(); i++) {
0808 const std::string& hltName = triggerNames.triggerName(i);
0809
0810
0811 for (unsigned int ipath = 0; ipath < triggerPath_.size(); ipath++) {
0812 if (hltName.find(triggerPath_[ipath]) != std::string::npos) {
0813 unsigned int triggerIndex(hltConfig.triggerIndex(hltName));
0814
0815
0816 if (triggerIndex < triggerResults->size()) {
0817 isFired = triggerResults->accept(triggerIndex);
0818 if (debug_ > 0)
0819 std::cout << triggerPath_[ipath] << " " << hltName << " " << isFired << std::endl;
0820 }
0821 }
0822 }
0823 }
0824 }
0825
0826 if (negateTrigger_ && isFired)
0827 return kContinue;
0828 else if (!(negateTrigger_) && !isFired)
0829 return kContinue;
0830
0831 #ifdef USE_CALLGRIND
0832 CALLGRIND_START_INSTRUMENTATION;
0833 #endif
0834
0835 if (debug_ > 0) {
0836 std::cout << "[MuScleFit-duringLoop]: loopCounter = " << loopCounter << " Run: " << event.id().run()
0837 << " Event: " << event.id().event() << std::endl;
0838 }
0839
0840
0841
0842
0843
0844
0845 if (loopCounter == 0) {
0846 if (!fastLoop || inputRootTreeFileName_.empty()) {
0847 if (debug_ > 0)
0848 std::cout << "Reading from edm event" << std::endl;
0849 selectMuons(event);
0850 duringFastLoop();
0851 ++totalEvents_;
0852 }
0853 }
0854
0855 return kContinue;
0856
0857 #ifdef USE_CALLGRIND
0858 CALLGRIND_STOP_INSTRUMENTATION;
0859 CALLGRIND_DUMP_STATS;
0860 #endif
0861 }
0862
0863 void MuScleFit::selectMuons(const edm::Event& event) {
0864 recMu1 = reco::Particle::LorentzVector(0, 0, 0, 0);
0865 recMu2 = reco::Particle::LorentzVector(0, 0, 0, 0);
0866
0867 std::vector<MuScleFitMuon> muons;
0868 muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
0869
0870
0871 if (debug_ > 0) {
0872 std::cout << "[MuScleFit::selectMuons] Debugging muons collections after call to muonSelector_->selectMuons"
0873 << std::endl;
0874 int iMu = 0;
0875 for (std::vector<MuScleFitMuon>::const_iterator it = muons.begin(); it < muons.end(); ++it) {
0876 std::cout << " - muon n. " << iMu << " = " << (*it) << std::endl;
0877 ++iMu;
0878 }
0879 }
0880
0881
0882
0883 std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes = MuScleFitUtils::findBestRecoRes(muons);
0884
0885 if (MuScleFitUtils::ResFound) {
0886 if (debug_ > 0) {
0887 std::cout << std::setprecision(9) << "Pt after findbestrecores: " << (recMuFromBestRes.first).Pt() << " "
0888 << (recMuFromBestRes.second).Pt() << std::endl;
0889 std::cout << "recMu1 = " << recMu1 << std::endl;
0890 std::cout << "recMu2 = " << recMu2 << std::endl;
0891 }
0892 recMu1 = recMuFromBestRes.first.p4();
0893 recMu2 = recMuFromBestRes.second.p4();
0894 recMuScleMu1 = recMuFromBestRes.first;
0895 recMuScleMu2 = recMuFromBestRes.second;
0896
0897 if (debug_ > 0) {
0898 std::cout << "after recMu1 = " << recMu1 << std::endl;
0899 std::cout << "after recMu2 = " << recMu2 << std::endl;
0900 std::cout << "mu1.pt = " << recMu1.Pt() << std::endl;
0901 std::cout << "mu2.pt = " << recMu2.Pt() << std::endl;
0902 std::cout << "after recMuScleMu1 = " << recMuScleMu1 << std::endl;
0903 std::cout << "after recMuScleMu2 = " << recMuScleMu2 << std::endl;
0904 }
0905 MuScleFitUtils::SavedPair.push_back(std::make_pair(recMu1, recMu2));
0906 MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(recMuScleMu1, recMuScleMu2));
0907 } else {
0908 MuScleFitUtils::SavedPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
0909 MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(MuScleFitMuon(), MuScleFitMuon()));
0910 }
0911
0912
0913
0914 UInt_t the_NVtx(0);
0915 Int_t the_numPUvtx(0);
0916 Float_t the_TrueNumInteractions(0);
0917
0918
0919
0920 edm::Handle<std::vector<PileupSummaryInfo> > puInfo;
0921 event.getByLabel(puInfoSrc_, puInfo);
0922 if (puInfo.isValid()) {
0923 std::vector<PileupSummaryInfo>::const_iterator PVI;
0924 for (PVI = puInfo->begin(); PVI != puInfo->end(); ++PVI) {
0925 int BX = PVI->getBunchCrossing();
0926 if (BX == 0) {
0927 the_TrueNumInteractions = PVI->getTrueNumInteractions();
0928 the_numPUvtx = PVI->getPU_NumInteractions();
0929 }
0930 }
0931 }
0932
0933 edm::Handle<std::vector<reco::Vertex> > vertices;
0934 event.getByLabel(vertexSrc_, vertices);
0935 if (vertices.isValid()) {
0936 std::vector<reco::Vertex>::const_iterator itv;
0937
0938 for (itv = vertices->begin(); itv != vertices->end(); ++itv) {
0939
0940 if (itv->ndof() < 5)
0941 continue;
0942 if (fabs(itv->z()) > 50.0)
0943 continue;
0944 if (fabs(itv->position().rho()) > 2.0)
0945 continue;
0946 ++the_NVtx;
0947 }
0948 }
0949
0950
0951 edm::Handle<GenEventInfoProduct> genEvtInfo;
0952 event.getByLabel("generator", genEvtInfo);
0953 double the_genEvtweight = 1.;
0954 if (genEvtInfo.isValid()) {
0955 the_genEvtweight = genEvtInfo->weight();
0956 }
0957
0958 muonPairs_.push_back(MuonPair(
0959 MuScleFitUtils::SavedPairMuScleFitMuons.back().first,
0960 MuScleFitUtils::SavedPairMuScleFitMuons.back().second,
0961 MuScleFitEvent(
0962 event.run(), event.id().event(), the_genEvtweight, the_numPUvtx, the_TrueNumInteractions, the_NVtx)));
0963
0964 if (MuScleFitUtils::speedup == false) {
0965 MuScleFitUtils::genPair.push_back(std::make_pair(genMuonPairs_.back().mu1.p4(), genMuonPairs_.back().mu2.p4()));
0966 }
0967 }
0968
0969 void MuScleFit::selectMuons(const int maxEvents, const TString& treeFileName) {
0970 std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
0971 RootTreeHandler rootTreeHandler;
0972 std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0973 if (MuScleFitUtils::speedup) {
0974 rootTreeHandler.readTree(
0975 maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPairMuScleFitMuons), theMuonType_, &evtRun);
0976 } else {
0977 rootTreeHandler.readTree(maxEvents,
0978 inputRootTreeFileName_,
0979 &(MuScleFitUtils::SavedPairMuScleFitMuons),
0980 theMuonType_,
0981 &evtRun,
0982 &(MuScleFitUtils::genMuscleFitPair));
0983 }
0984
0985 std::vector<std::pair<unsigned int, unsigned long long> >::iterator evtRunIt = evtRun.begin();
0986 std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator it = MuScleFitUtils::SavedPairMuScleFitMuons.begin();
0987 std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator genIt;
0988 if (MuScleFitUtils::speedup == false)
0989 genIt = MuScleFitUtils::genMuscleFitPair.begin();
0990 for (; it != MuScleFitUtils::SavedPairMuScleFitMuons.end(); ++it, ++evtRunIt) {
0991
0992
0993
0994 double pt1 = it->first.pt();
0995
0996 double pt2 = it->second.pt();
0997
0998 double eta1 = it->first.eta();
0999
1000 double eta2 = it->second.eta();
1001
1002
1003 bool dontPass = false;
1004 bool eta1InFirstRange;
1005 bool eta2InFirstRange;
1006 bool eta1InSecondRange;
1007 bool eta2InSecondRange;
1008
1009 int ch1 = it->first.charge();
1010 int ch2 = it->second.charge();
1011
1012 if (MuScleFitUtils::separateRanges_) {
1013 eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1014 eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1015 eta1InSecondRange =
1016 eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1017 eta2InSecondRange =
1018 eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1019
1020
1021 if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1022 pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1023 ((eta1InFirstRange && eta2InSecondRange && ch1 >= ch2) ||
1024 (eta1InSecondRange && eta2InFirstRange && ch1 < ch2))))
1025 dontPass = true;
1026 } else {
1027 eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1028 eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1029 eta1InSecondRange =
1030 eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1031 eta2InSecondRange =
1032 eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1033 if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1034 pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1035 (((eta1InFirstRange && !eta2InFirstRange) && (eta2InSecondRange && !eta1InSecondRange) && ch1 >= ch2) ||
1036 ((eta2InFirstRange && !eta1InFirstRange) && (eta1InSecondRange && !eta2InSecondRange) && ch1 < ch2))))
1037 dontPass = true;
1038 }
1039
1040
1041 double deltaPhi = MuScleFitUtils::deltaPhi(it->first.phi(), it->second.phi());
1042 if ((deltaPhi <= MuScleFitUtils::deltaPhiMinCut_) || (deltaPhi >= MuScleFitUtils::deltaPhiMaxCut_))
1043 dontPass = true;
1044
1045 lorentzVector vec1 = it->first.p4();
1046 lorentzVector vec2 = it->second.p4();
1047 if (ch1 >= ch2) {
1048 lorentzVector vectemp = vec1;
1049 vec1 = vec2;
1050 vec2 = vectemp;
1051 }
1052
1053 if (!dontPass) {
1054
1055 if ((MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0)) {
1056 applySmearing(vec1);
1057 applyBias(vec1, -1);
1058 applySmearing(vec2);
1059 applyBias(vec2, 1);
1060 }
1061
1062 MuScleFitUtils::SavedPair.push_back(std::make_pair(vec1, vec2));
1063 }
1064
1065
1066 muonPairs_.push_back(MuonPair(
1067 MuScleFitMuon(vec1, -1),
1068 MuScleFitMuon(vec2, +1),
1069 MuScleFitEvent(
1070 (*evtRunIt).first, (*evtRunIt).second, 0, 0, 0, 0))
1071 );
1072
1073
1074 if (!MuScleFitUtils::speedup) {
1075 MuScleFitUtils::genPair.push_back(std::make_pair(genIt->first.p4(), genIt->second.p4()));
1076 genMuonPairs_.push_back(GenMuonPair(genIt->first.p4(), genIt->second.p4(), 0));
1077 ++genIt;
1078 }
1079 }
1080 plotter->fillTreeRec(MuScleFitUtils::SavedPair);
1081 if (!(MuScleFitUtils::speedup)) {
1082 plotter->fillTreeGen(MuScleFitUtils::genPair);
1083 }
1084 }
1085
1086 void MuScleFit::duringFastLoop() {
1087
1088
1089 MuScleFitUtils::ResFound = false;
1090 recMu1 = (MuScleFitUtils::SavedPair[iev].first);
1091 recMu2 = (MuScleFitUtils::SavedPair[iev].second);
1092
1093
1094
1095 if (recMu1.Pt() > 0 && recMu2.Pt() > 0) {
1096 MuScleFitUtils::ResFound = true;
1097 if (debug_ > 0)
1098 std::cout << "Ev = " << iev << ": found muons in tree with Pt = " << recMu1.Pt() << " " << recMu2.Pt()
1099 << std::endl;
1100 }
1101
1102 if (debug_ > 0)
1103 std::cout << "About to start lik par correction and histo filling; ResFound is " << MuScleFitUtils::ResFound
1104 << std::endl;
1105
1106
1107 if (MuScleFitUtils::ResFound) {
1108
1109
1110
1111
1112 double weight = MuScleFitUtils::computeWeight((recMu1 + recMu2).mass(), iev, true);
1113 if (debug_ > 0) {
1114 std::cout << "Loop #" << loopCounter << "Event #" << iev << ": before correction Pt1 = " << recMu1.Pt()
1115 << " Pt2 = " << recMu2.Pt() << std::endl;
1116 }
1117
1118
1119 if (loopCounter > 0) {
1120 if (MuScleFitUtils::doScaleFit[loopCounter - 1]) {
1121 recMu1 = (MuScleFitUtils::applyScale(recMu1, MuScleFitUtils::parvalue[loopCounter - 1], -1));
1122 recMu2 = (MuScleFitUtils::applyScale(recMu2, MuScleFitUtils::parvalue[loopCounter - 1], 1));
1123 }
1124 }
1125 if (debug_ > 0) {
1126 std::cout << "Loop #" << loopCounter << "Event #" << iev << ": after correction Pt1 = " << recMu1.Pt()
1127 << " Pt2 = " << recMu2.Pt() << std::endl;
1128 }
1129
1130 reco::Particle::LorentzVector bestRecRes(recMu1 + recMu2);
1131
1132
1133
1134
1135 mapHisto_["hRecBestMu"]->Fill(recMu1, -1, weight);
1136 mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
1137 mapHisto_["hRecBestMu"]->Fill(recMu2, +1, weight);
1138 mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
1139 mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
1140
1141 mapHisto_["hRecBestRes"]->Fill(bestRecRes, +1, weight);
1142 mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, +1, 1.);
1143
1144
1145
1146
1147
1148
1149 mapHisto_["hRecBestResVSMu"]->Fill(recMu1, bestRecRes, -1, weight);
1150 mapHisto_["hRecBestResVSMu"]->Fill(recMu2, bestRecRes, +1, weight);
1151
1152 mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes, weight);
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164 mapHisto_["hRecBestResVSRes"]->Fill(bestRecRes, bestRecRes, +1, weight);
1165
1166 std::vector<double>* parval;
1167 std::vector<double> initpar;
1168
1169
1170 if (loopCounter == 0) {
1171 initpar = MuScleFitUtils::parResol;
1172 initpar.insert(initpar.end(), MuScleFitUtils::parScale.begin(), MuScleFitUtils::parScale.end());
1173 initpar.insert(initpar.end(), MuScleFitUtils::parCrossSection.begin(), MuScleFitUtils::parCrossSection.end());
1174 initpar.insert(initpar.end(), MuScleFitUtils::parBgr.begin(), MuScleFitUtils::parBgr.end());
1175 parval = &initpar;
1176 } else {
1177 parval = &(MuScleFitUtils::parvalue[loopCounter - 1]);
1178 }
1179
1180
1181
1182 if (!MuScleFitUtils::speedup) {
1183
1184 if (checkDeltaR(MuScleFitUtils::genPair[iev].first, recMu1)) {
1185 fillComparisonHistograms(MuScleFitUtils::genPair[iev].first, recMu1, "Gen", -1);
1186 }
1187 if (checkDeltaR(MuScleFitUtils::genPair[iev].second, recMu2)) {
1188 fillComparisonHistograms(MuScleFitUtils::genPair[iev].second, recMu2, "Gen", +1);
1189 }
1190 if (compareToSimTracks_) {
1191
1192 if (checkDeltaR(MuScleFitUtils::simPair[iev].first, recMu1)) {
1193 fillComparisonHistograms(MuScleFitUtils::simPair[iev].first, recMu1, "Sim", -1);
1194 }
1195 if (checkDeltaR(MuScleFitUtils::simPair[iev].second, recMu2)) {
1196 fillComparisonHistograms(MuScleFitUtils::simPair[iev].second, recMu2, "Sim", +1);
1197 }
1198 }
1199 }
1200
1201
1202
1203
1204
1205 mapHisto_["hFunctionResolPt"]->Fill(
1206 recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1207 mapHisto_["hFunctionResolCotgTheta"]->Fill(
1208 recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1209 mapHisto_["hFunctionResolPhi"]->Fill(
1210 recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1211 mapHisto_["hFunctionResolPt"]->Fill(
1212 recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1213 mapHisto_["hFunctionResolCotgTheta"]->Fill(
1214 recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1215 mapHisto_["hFunctionResolPhi"]->Fill(
1216 recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1217
1218
1219
1220 if (debug_ > 0)
1221 std::cout << "mass = " << bestRecRes.mass() << std::endl;
1222 if (weight != 0.) {
1223 double massResol;
1224 double prob;
1225 double deltalike;
1226 if (loopCounter == 0) {
1227 std::vector<double> initpar;
1228 initpar.reserve((int)(MuScleFitUtils::parResol.size()));
1229 for (int i = 0; i < (int)(MuScleFitUtils::parResol.size()); i++) {
1230 initpar.push_back(MuScleFitUtils::parResol[i]);
1231 }
1232 for (int i = 0; i < (int)(MuScleFitUtils::parScale.size()); i++) {
1233 initpar.push_back(MuScleFitUtils::parScale[i]);
1234 }
1235
1236
1237
1238 MuScleFitUtils::crossSectionHandler->addParameters(initpar);
1239
1240 for (int i = 0; i < (int)(MuScleFitUtils::parBgr.size()); i++) {
1241 initpar.push_back(MuScleFitUtils::parBgr[i]);
1242 }
1243 massResol = MuScleFitUtils::massResolution(recMu1, recMu2, initpar);
1244
1245 prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1246 bestRecRes.Eta(),
1247 bestRecRes.Rapidity(),
1248 massResol,
1249 initpar,
1250 true,
1251 recMu1.eta(),
1252 recMu2.eta());
1253 } else {
1254 massResol = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parvalue[loopCounter - 1]);
1255
1256
1257 prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1258 bestRecRes.Eta(),
1259 bestRecRes.Rapidity(),
1260 massResol,
1261 MuScleFitUtils::parvalue[loopCounter - 1],
1262 true,
1263 recMu1.eta(),
1264 recMu2.eta());
1265 }
1266 if (debug_ > 0)
1267 std::cout << "inside weight: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1268 if (prob > 0) {
1269 if (debug_ > 0)
1270 std::cout << "inside prob: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1271
1272 deltalike = log(prob) * weight;
1273 mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike);
1274 mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike);
1275 mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike);
1276 mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike);
1277
1278 double recoMass = (recMu1 + recMu2).mass();
1279 if (recoMass != 0) {
1280
1281 mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1);
1282 mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1);
1283 mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol / recoMass, -1);
1284 mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol / recoMass, +1);
1285 }
1286
1287 if (MuScleFitUtils::debugMassResol_) {
1288 mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1);
1289 mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1);
1290 mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1);
1291 mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1);
1292 mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1);
1293 mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1);
1294 }
1295
1296 if (!MuScleFitUtils::speedup) {
1297 double genMass = (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second).mass();
1298
1299 if (genMass != 0) {
1300 mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].first),
1301 (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1302 -1);
1303 mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].second),
1304 (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1305 +1);
1306 double diffMass = (recoMass - genMass) / genMass;
1307
1308
1309 double pt1 = recMu1.pt();
1310 double eta1 = recMu1.eta();
1311 double pt2 = recMu2.pt();
1312 double eta2 = recMu2.eta();
1313
1314 if (diffMass == diffMass) {
1315
1316 mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
1317 mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
1318 mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
1319 mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
1320
1321 mapHisto_["hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
1322 mapHisto_["hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
1323 } else {
1324 std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
1325 }
1326 }
1327
1328 double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
1329 mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
1330 mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
1331 }
1332
1333 mapHisto_["hMass_P"]->Fill(bestRecRes.mass(), prob);
1334 if (debug_ > 0)
1335 std::cout << "mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1336 mapHisto_["hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
1337
1338 mapHisto_["hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
1339 mapHisto_["hMassProbVsMu"]->Fill(recMu1, bestRecRes, -1, prob);
1340 mapHisto_["hMassProbVsMu"]->Fill(recMu2, bestRecRes, +1, prob);
1341 mapHisto_["hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
1342 mapHisto_["hMassProbVsMu_fine"]->Fill(recMu1, bestRecRes, -1, prob);
1343 mapHisto_["hMassProbVsMu_fine"]->Fill(recMu2, bestRecRes, +1, prob);
1344 }
1345 }
1346 }
1347
1348
1349
1350 if (loopCounter > 0) {
1351 if (debug_ > 0)
1352 std::cout << "[MuScleFit]: filling the pair" << std::endl;
1353 MuScleFitUtils::SavedPair[iev] = std::make_pair(recMu1, recMu2);
1354 }
1355
1356 iev++;
1357 MuScleFitUtils::iev_++;
1358
1359
1360 }
1361
1362 bool MuScleFit::checkDeltaR(reco::Particle::LorentzVector& genMu, reco::Particle::LorentzVector& recMu) {
1363
1364 double deltaR =
1365 sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
1366 ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
1367 if (deltaR < 0.01)
1368 return true;
1369 else if (debug_ > 0) {
1370 std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
1371 << " DOES NOT MATCH with generated muon from resonance: " << std::endl
1372 << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
1373 }
1374 return false;
1375 }
1376
1377 void MuScleFit::fillComparisonHistograms(const reco::Particle::LorentzVector& genMu,
1378 const reco::Particle::LorentzVector& recMu,
1379 const std::string& inputName,
1380 const int charge) {
1381 std::string name(inputName + "VSMu");
1382 mapHisto_["hResolPt" + name]->Fill(recMu, (-genMu.Pt() + recMu.Pt()) / genMu.Pt(), charge);
1383 mapHisto_["hResolTheta" + name]->Fill(recMu, (-genMu.Theta() + recMu.Theta()), charge);
1384 mapHisto_["hResolCotgTheta" + name]->Fill(
1385 recMu, (-cos(genMu.Theta()) / sin(genMu.Theta()) + cos(recMu.Theta()) / sin(recMu.Theta())), charge);
1386 mapHisto_["hResolEta" + name]->Fill(recMu, (-genMu.Eta() + recMu.Eta()), charge);
1387 mapHisto_["hResolPhi" + name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);
1388
1389
1390 if ((genMu.Pt() != 0) && (recMu.Pt() != 0)) {
1391 mapHisto_["hPtRecoVsPt" + inputName]->Fill(genMu.Pt(), recMu.Pt());
1392 }
1393 }
1394
1395 void MuScleFit::applySmearing(reco::Particle::LorentzVector& mu) {
1396 if (MuScleFitUtils::SmearType > 0) {
1397 mu = MuScleFitUtils::applySmearing(mu);
1398 if (debug_ > 0)
1399 std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after smearing Pt = " << mu.Pt() << std::endl;
1400 }
1401 }
1402
1403 void MuScleFit::applyBias(reco::Particle::LorentzVector& mu, const int charge) {
1404 if (MuScleFitUtils::BiasType > 0) {
1405 mu = MuScleFitUtils::applyBias(mu, charge);
1406 if (debug_ > 0)
1407 std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after bias Pt = " << mu.Pt() << std::endl;
1408 }
1409 }
1410
1411
1412
1413 void MuScleFit::checkParameters() {
1414
1415 if (MuScleFitUtils::doResolFit.size() != maxLoopNumber) {
1416 std::cout << "[MuScleFit-Constructor]: wrong size of resolution fits selector = "
1417 << MuScleFitUtils::doResolFit.size() << std::endl;
1418 std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1419 abort();
1420 }
1421 if (MuScleFitUtils::doScaleFit.size() != maxLoopNumber) {
1422 std::cout << "[MuScleFit-Constructor]: wrong size of scale fits selector = " << MuScleFitUtils::doScaleFit.size()
1423 << std::endl;
1424 std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1425 abort();
1426 }
1427 if (MuScleFitUtils::doCrossSectionFit.size() != maxLoopNumber) {
1428 std::cout << "[MuScleFit-Constructor]: wrong size of cross section fits selector = "
1429 << MuScleFitUtils::doCrossSectionFit.size() << std::endl;
1430 std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1431 abort();
1432 }
1433 if (MuScleFitUtils::doBackgroundFit.size() != maxLoopNumber) {
1434 std::cout << "[MuScleFit-Constructor]: wrong size of background fits selector = "
1435 << MuScleFitUtils::doBackgroundFit.size() << std::endl;
1436 std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1437 abort();
1438 }
1439
1440
1441
1442 if ((MuScleFitUtils::BiasType == 1 && MuScleFitUtils::parBias.size() != 2) ||
1443 (MuScleFitUtils::BiasType == 2 && MuScleFitUtils::parBias.size() != 2) ||
1444 (MuScleFitUtils::BiasType == 3 && MuScleFitUtils::parBias.size() != 4) ||
1445 (MuScleFitUtils::BiasType == 4 && MuScleFitUtils::parBias.size() != 3) ||
1446 (MuScleFitUtils::BiasType == 5 && MuScleFitUtils::parBias.size() != 3) ||
1447 (MuScleFitUtils::BiasType == 6 &&
1448 MuScleFitUtils::parBias.size() != 3) ||
1449 (MuScleFitUtils::BiasType == 7 &&
1450 MuScleFitUtils::parBias.size() != 4) ||
1451 (MuScleFitUtils::BiasType == 8 && MuScleFitUtils::parBias.size() != 4) ||
1452 (MuScleFitUtils::BiasType == 9 && MuScleFitUtils::parBias.size() != 2) ||
1453 (MuScleFitUtils::BiasType == 10 && MuScleFitUtils::parBias.size() != 3) ||
1454 (MuScleFitUtils::BiasType == 11 &&
1455 MuScleFitUtils::parBias.size() != 4) ||
1456 (MuScleFitUtils::BiasType == 12 &&
1457 MuScleFitUtils::parBias.size() != 6) ||
1458 (MuScleFitUtils::BiasType == 13 &&
1459 MuScleFitUtils::parBias.size() != 8) ||
1460 MuScleFitUtils::BiasType < 0 ||
1461 MuScleFitUtils::BiasType > 13) {
1462 std::cout << "[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
1463 abort();
1464 }
1465
1466
1467 if ((MuScleFitUtils::SmearType == 1 && MuScleFitUtils::parSmear.size() != 3) ||
1468 (MuScleFitUtils::SmearType == 2 && MuScleFitUtils::parSmear.size() != 4) ||
1469 (MuScleFitUtils::SmearType == 3 && MuScleFitUtils::parSmear.size() != 5) ||
1470 (MuScleFitUtils::SmearType == 4 && MuScleFitUtils::parSmear.size() != 6) ||
1471 (MuScleFitUtils::SmearType == 5 && MuScleFitUtils::parSmear.size() != 7) ||
1472 (MuScleFitUtils::SmearType == 6 && MuScleFitUtils::parSmear.size() != 16) ||
1473 (MuScleFitUtils::SmearType == 7 && !MuScleFitUtils::parSmear.empty()) || MuScleFitUtils::SmearType < 0 ||
1474 MuScleFitUtils::SmearType > 7) {
1475 std::cout << "[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
1476 abort();
1477 }
1478
1479
1480 if (MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolFix.size() ||
1481 MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolOrder.size()) {
1482 std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
1483 abort();
1484 }
1485 if (MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleFix.size() ||
1486 MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleOrder.size()) {
1487 std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
1488 abort();
1489 }
1490 if (MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionFix.size() ||
1491 MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionOrder.size()) {
1492 std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1493 abort();
1494 }
1495 if (MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrFix.size() ||
1496 MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrOrder.size()) {
1497 std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1498 abort();
1499 }
1500
1501
1502
1503 if (MuScleFitUtils::resfind.size() != 6) {
1504 std::cout << "[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
1505 abort();
1506 }
1507 }
1508
1509 bool MuScleFit::selGlobalMuon(const pat::Muon* aMuon) {
1510 reco::TrackRef iTrack = aMuon->innerTrack();
1511 const reco::HitPattern& p = iTrack->hitPattern();
1512
1513 reco::TrackRef gTrack = aMuon->globalTrack();
1514 const reco::HitPattern& q = gTrack->hitPattern();
1515
1516 return (
1517 iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 && q.numberOfValidMuonHits() > 0 &&
1518 iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
1519 aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
1520 fabs(iTrack->dxy()) < 3.0 &&
1521 fabs(iTrack->dz()) < 15.0);
1522 }
1523
1524 bool MuScleFit::selTrackerMuon(const pat::Muon* aMuon) {
1525 reco::TrackRef iTrack = aMuon->innerTrack();
1526 const reco::HitPattern& p = iTrack->hitPattern();
1527
1528 return (
1529 iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
1530 aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
1531 fabs(iTrack->dxy()) < 3.0 &&
1532 fabs(iTrack->dz()) < 15.0);
1533 }
1534
1535 #include "FWCore/Framework/interface/MakerMacros.h"
1536 DEFINE_FWK_LOOPER(MuScleFit);