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