Inheritance diagram for MuScleFit:

Public Member Functions
void	beginOfJobInConstructor ()

virtual void	duringFastLoop ()

edm::EDLooper::Status	duringLoop (const edm::Event &event, const edm::EventSetup &eventSetup) override

virtual void	endOfFastLoop (const unsigned int iLoop)

void	endOfJob () override

edm::EDLooper::Status	endOfLoop (const edm::EventSetup &eventSetup, unsigned int iLoop) override

template<typename T >
std::vector< MuScleFitMuon >	fillMuonCollection (const std::vector< T > &tracks)

	MuScleFit (const edm::ParameterSet &pset)

void	startingNewLoop (unsigned int iLoop) override

	~MuScleFit () override

Public Member Functions inherited from edm::EDLooper
	EDLooper ()

	EDLooper (EDLooper const &)=delete

EDLooper &	operator= (EDLooper const &)=delete

	~EDLooper () override

Public Member Functions inherited from edm::EDLooperBase
virtual void	attachTo (ActivityRegistry &)
	Override this method if you need to monitor the state of the processing. More...

void	beginOfJob (EventSetupImpl const &)

virtual void	beginOfJob (EventSetup const &)

virtual void	beginOfJob ()

void	copyInfo (ScheduleInfo const &)

void	doBeginLuminosityBlock (LuminosityBlockPrincipal &, EventSetupImpl const &, ProcessContext *)

void	doBeginRun (RunPrincipal &, EventSetupImpl const &, ProcessContext *)

Status	doDuringLoop (EventPrincipal &eventPrincipal, EventSetupImpl const &es, ProcessingController &, StreamContext *)

void	doEndLuminosityBlock (LuminosityBlockPrincipal &, EventSetupImpl const &, ProcessContext *)

Status	doEndOfLoop (EventSetupImpl const &es)

void	doEndRun (RunPrincipal &, EventSetupImpl const &, ProcessContext *)

void	doStartingNewLoop ()

	EDLooperBase ()

	EDLooperBase (EDLooperBase const &)=delete

virtual std::set< eventsetup::EventSetupRecordKey >	modifyingRecords () const

EDLooperBase &	operator= (EDLooperBase const &)=delete

void	prepareForNextLoop (eventsetup::EventSetupProvider *esp)

void	setActionTable (ExceptionToActionTable const *actionTable)

void	setModuleChanger (ModuleChanger *)

virtual	~EDLooperBase () noexcept(false)

Protected Member Functions
void	applyBias (reco::Particle::LorentzVector &mu, const int charge)
	Apply the bias if needed using the function in MuScleFitUtils. More...

void	applySmearing (reco::Particle::LorentzVector &mu)
	Apply the smearing if needed using the function in MuScleFitUtils. More...

bool	checkDeltaR (reco::Particle::LorentzVector &genMu, reco::Particle::LorentzVector &recMu)
	Check if two lorentzVector are near in deltaR. More...

void	checkParameters ()

void	fillComparisonHistograms (const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recoMu, const std::string &inputName, const int charge)
	Fill the reco vs gen and reco vs sim comparison histograms. More...

void	selectMuons (const edm::Event &event)

void	selectMuons (const int maxEvents, const TString &treeFileName)

bool	selGlobalMuon (const pat::Muon *aMuon)
	Function for onia selections. More...

bool	selTrackerMuon (const pat::Muon *aMuon)

template<typename T >
void	takeSelectedMuonType (const T &muon, std::vector< reco::Track > &tracks)
	Template method used to fill the track collection starting from reco::muons or pat::muons. More...

Protected Member Functions inherited from edm::EDLooperBase
ModuleChanger *	moduleChanger ()
	This only returns a non-zero value during the call to endOfLoop. More...

ScheduleInfo const *	scheduleInfo () const
	This returns a non-zero value after the constructor has been called. More...

Protected Attributes
bool	compareToSimTracks_

bool	fastLoop

std::string	genParticlesName_

int	iev

bool	ifGenPart

bool	ifHepMC

std::string	inputRootTreeFileName_

unsigned int	loopCounter

int	maxEventsFromRootTree_

unsigned int	maxLoopNumber

double	maxResMass_hwindow [6]

double	minResMass_hwindow [6]

std::unique_ptr< MuScleFitMuonSelector >	muonSelector_

bool	negateTrigger_

int	numberOfEwkZ

int	numberOfSimMuons

int	numberOfSimTracks

int	numberOfSimVertices

std::string	outputRootTreeFileName_

bool	PATmuons_

MuScleFitPlotter *	plotter

edm::InputTag	puInfoSrc_

reco::Particle::LorentzVector	recMu1

reco::Particle::LorentzVector	recMu2

MuScleFitMuon	recMuScleMu1

MuScleFitMuon	recMuScleMu2

bool	saveAllToTree_

edm::InputTag	simTracksCollection_

MuonServiceProxy *	theService

int	totalEvents_

std::vector< std::string >	triggerPath_

std::string	triggerResultsLabel_

std::string	triggerResultsProcess_

edm::InputTag	vertexSrc_

Additional Inherited Members
Public Types inherited from edm::EDLooperBase
enum	Status { kContinue, kStop }

Private Member Functions inherited from MuScleFitBase
	MuScleFitBase (const edm::ParameterSet &iConfig)

virtual	~MuScleFitBase ()(false)

void	clearHistoMap ()
	Clean the histograms map. More...

void	fillHistoMap (TFile *outputFile, unsigned int iLoop)
	Create the histograms map. More...

void	readProbabilityDistributionsFromFile ()
	Read probability distributions from a local root file. More...

void	writeHistoMap (const unsigned int iLoop)
	Save the histograms map to file. More...

Private Attributes inherited from MuScleFitBase
int	debug_

std::vector< GenMuonPair >	genMuonPairs_
	Stores the genMuon pairs and the motherId prior to the creation of the internal tree. More...

std::map< std::string, Histograms * >	mapHisto_
	The map of histograms. More...

std::vector< MuonPair >	muonPairs_
	Used to store the muon pairs plus run and event number prior to the creation of the internal tree. More...

std::string	probabilitiesFile_

std::string	probabilitiesFileInPath_

std::vector< TFile * >	theFiles_
	The files were the histograms are saved. More...

std::string	theGenInfoRootFileName_

edm::InputTag	theMuonLabel_

int	theMuonType_

std::string	theRootFileName_

Detailed Description

Analyzer of the Global muon tracks

Definition at line 183 of file MuScleFit.cc.

Constructor & Destructor Documentation

MuScleFit::MuScleFit ( const edm::ParameterSet & pset )

Definition at line 371 of file MuScleFit.cc.

                                                 : MuScleFitBase(pset), totalEvents_(0) {
   MuScleFitUtils::debug = debug_;
   if (debug_ > 0)
     std::cout << "[MuScleFit]: Constructor" << std::endl;
 
   if ((theMuonType_ < -4 || theMuonType_ > 5) && theMuonType_ < 10) {
     std::cout << "[MuScleFit]: Unknown muon type! Aborting." << std::endl;
     abort();
   }
 
   loopCounter = 0;
 
   // Boundaries for h-function computation (to be improved!)
   // -------------------------------------------------------
   minResMass_hwindow[0] = 71.1876;  // 76.;
   maxResMass_hwindow[0] = 111.188;  // 106.;
   minResMass_hwindow[1] = 10.15;
   maxResMass_hwindow[1] = 10.55;
   minResMass_hwindow[2] = 9.8;
   maxResMass_hwindow[2] = 10.2;
   minResMass_hwindow[3] = 9.25;
   maxResMass_hwindow[3] = 9.65;
   minResMass_hwindow[4] = 3.58;
   maxResMass_hwindow[4] = 3.78;
   minResMass_hwindow[5] = 3.0;
   maxResMass_hwindow[5] = 3.2;
 
   // Max number of loops (if > 2 then try to minimize likelihood more than once)
   // ---------------------------------------------------------------------------
   maxLoopNumber = pset.getUntrackedParameter<int>("maxLoopNumber", 2);
   fastLoop = pset.getUntrackedParameter<bool>("FastLoop", true);
 
   // Selection of fits according to loop
   MuScleFitUtils::doResolFit = pset.getParameter<std::vector<int> >("doResolFit");
   MuScleFitUtils::doScaleFit = pset.getParameter<std::vector<int> >("doScaleFit");
   MuScleFitUtils::doCrossSectionFit = pset.getParameter<std::vector<int> >("doCrossSectionFit");
   MuScleFitUtils::doBackgroundFit = pset.getParameter<std::vector<int> >("doBackgroundFit");
 
   // Bias and smear types
   // --------------------
   int biasType = pset.getParameter<int>("BiasType");
   MuScleFitUtils::BiasType = biasType;
   // No error, the scale functions are used also for the bias
   MuScleFitUtils::biasFunction = scaleFunctionVecService(biasType);
   int smearType = pset.getParameter<int>("SmearType");
   MuScleFitUtils::SmearType = smearType;
   MuScleFitUtils::smearFunction = smearFunctionService(smearType);
 
   // Fit types
   // ---------
   int resolFitType = pset.getParameter<int>("ResolFitType");
   MuScleFitUtils::ResolFitType = resolFitType;
   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType);
   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolFitType);
   int scaleType = pset.getParameter<int>("ScaleFitType");
   MuScleFitUtils::ScaleFitType = scaleType;
   MuScleFitUtils::scaleFunction = scaleFunctionService(scaleType);
   MuScleFitUtils::scaleFunctionForVec = scaleFunctionVecService(scaleType);
 
   // Initial parameters values
   // -------------------------
   MuScleFitUtils::parBias = pset.getParameter<std::vector<double> >("parBias");
   MuScleFitUtils::parSmear = pset.getParameter<std::vector<double> >("parSmear");
   MuScleFitUtils::parResol = pset.getParameter<std::vector<double> >("parResol");
   MuScleFitUtils::parResolStep =
       pset.getUntrackedParameter<std::vector<double> >("parResolStep", std::vector<double>());
   MuScleFitUtils::parResolMin = pset.getUntrackedParameter<std::vector<double> >("parResolMin", std::vector<double>());
   MuScleFitUtils::parResolMax = pset.getUntrackedParameter<std::vector<double> >("parResolMax", std::vector<double>());
   MuScleFitUtils::parScale = pset.getParameter<std::vector<double> >("parScale");
   MuScleFitUtils::parScaleStep =
       pset.getUntrackedParameter<std::vector<double> >("parScaleStep", std::vector<double>());
   MuScleFitUtils::parScaleMin = pset.getUntrackedParameter<std::vector<double> >("parScaleMin", std::vector<double>());
   MuScleFitUtils::parScaleMax = pset.getUntrackedParameter<std::vector<double> >("parScaleMax", std::vector<double>());
   MuScleFitUtils::parCrossSection = pset.getParameter<std::vector<double> >("parCrossSection");
   MuScleFitUtils::parBgr = pset.getParameter<std::vector<double> >("parBgr");
   MuScleFitUtils::parResolFix = pset.getParameter<std::vector<int> >("parResolFix");
   MuScleFitUtils::parScaleFix = pset.getParameter<std::vector<int> >("parScaleFix");
   MuScleFitUtils::parCrossSectionFix = pset.getParameter<std::vector<int> >("parCrossSectionFix");
   MuScleFitUtils::parBgrFix = pset.getParameter<std::vector<int> >("parBgrFix");
   MuScleFitUtils::parResolOrder = pset.getParameter<std::vector<int> >("parResolOrder");
   MuScleFitUtils::parScaleOrder = pset.getParameter<std::vector<int> >("parScaleOrder");
   MuScleFitUtils::parCrossSectionOrder = pset.getParameter<std::vector<int> >("parCrossSectionOrder");
   MuScleFitUtils::parBgrOrder = pset.getParameter<std::vector<int> >("parBgrOrder");
 
   MuScleFitUtils::resfind = pset.getParameter<std::vector<int> >("resfind");
   MuScleFitUtils::FitStrategy = pset.getParameter<int>("FitStrategy");
 
   // Option to skip unnecessary stuff
   // --------------------------------
   MuScleFitUtils::speedup = pset.getParameter<bool>("speedup");
 
   // Option to skip simTracks comparison
   compareToSimTracks_ = pset.getParameter<bool>("compareToSimTracks");
   simTracksCollection_ = pset.getUntrackedParameter<edm::InputTag>("SimTracksCollection", edm::InputTag("g4SimHits"));
 
   triggerResultsLabel_ = pset.getUntrackedParameter<std::string>("TriggerResultsLabel");
   triggerResultsProcess_ = pset.getUntrackedParameter<std::string>("TriggerResultsProcess");
   triggerPath_ = pset.getUntrackedParameter<std::vector<std::string> >("TriggerPath");
   negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
   saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
 
   // input collections for PU related infos
   puInfoSrc_ = pset.getUntrackedParameter<edm::InputTag>("PileUpSummaryInfo");
   vertexSrc_ = pset.getUntrackedParameter<edm::InputTag>("PrimaryVertexCollection");
 
   PATmuons_ = pset.getUntrackedParameter<bool>("PATmuons", false);
   genParticlesName_ = pset.getUntrackedParameter<std::string>("GenParticlesName", "genParticles");
 
   // Use the probability file or not. If not it will perform a simpler selection taking the muon pair with
   // invariant mass closer to the pdf value and will crash if some fit is attempted.
   MuScleFitUtils::useProbsFile_ = pset.getUntrackedParameter<bool>("UseProbsFile", true);
 
   // This must be set to true if using events generated with Sherpa
   MuScleFitUtils::sherpa_ = pset.getUntrackedParameter<bool>("Sherpa", false);
 
   MuScleFitUtils::rapidityBinsForZ_ = pset.getUntrackedParameter<bool>("RapidityBinsForZ", true);
 
   // Set the cuts on muons to be used in the fit
   MuScleFitUtils::separateRanges_ = pset.getUntrackedParameter<bool>("SeparateRanges", true);
   MuScleFitUtils::maxMuonPt_ = pset.getUntrackedParameter<double>("MaxMuonPt", 100000000.);
   MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
   MuScleFitUtils::minMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MinMuonEtaFirstRange", -6.);
   MuScleFitUtils::maxMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaFirstRange", 6.);
   MuScleFitUtils::minMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MinMuonEtaSecondRange", -100.);
   MuScleFitUtils::maxMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaSecondRange", 100.);
   MuScleFitUtils::deltaPhiMinCut_ = pset.getUntrackedParameter<double>("DeltaPhiMinCut", -100.);
   MuScleFitUtils::deltaPhiMaxCut_ = pset.getUntrackedParameter<double>("DeltaPhiMaxCut", 100.);
 
   MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("DebugMassResol", false);
   // MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents;
 
   // Check for parameters consistency
   // it will abort in case of errors.
   checkParameters();
 
   // Generate array of gaussian-distributed numbers for smearing
   // -----------------------------------------------------------
   if (MuScleFitUtils::SmearType > 0) {
     std::cout << "[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
     TF1 G("G", "[0]*exp(-0.5*pow(x,2))", -5., 5.);
     double norm = 1 / sqrt(2 * TMath::Pi());
     G.SetParameter(0, norm);
     for (int i = 0; i < 10000; i++) {
       for (int j = 0; j < 7; j++) {
         MuScleFitUtils::x[j][i] = G.GetRandom();
       }
     }
   }
   MuScleFitUtils::goodmuon = 0;
 
   if (theMuonType_ > 0 && theMuonType_ < 4) {
     MuScleFitUtils::MuonTypeForCheckMassWindow = theMuonType_ - 1;
     MuScleFitUtils::MuonType = theMuonType_ - 1;
   } else if (theMuonType_ == 0 || theMuonType_ == 4 || theMuonType_ == 5 || theMuonType_ >= 10 || theMuonType_ == -1 ||
              theMuonType_ == -2 || theMuonType_ == -3 || theMuonType_ == -4) {
     MuScleFitUtils::MuonTypeForCheckMassWindow = 2;
     MuScleFitUtils::MuonType = 2;
   } else {
     std::cout << "Wrong muon type " << theMuonType_ << std::endl;
     exit(1);
   }
 
   // When using standalone muons switch to the single Z pdf
   if (theMuonType_ == 2) {
     MuScleFitUtils::rapidityBinsForZ_ = false;
   }
 
   // Initialize ResMaxSigma And ResHalfWidth - 0 = global, 1 = SM, 2 = tracker
   // -------------------------------------------------------------------------
   MuScleFitUtils::massWindowHalfWidth[0][0] = 20.;
   MuScleFitUtils::massWindowHalfWidth[0][1] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[0][2] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[0][3] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[0][4] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[0][5] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[1][0] = 20.;
   MuScleFitUtils::massWindowHalfWidth[1][1] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[1][2] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[1][3] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[1][4] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[1][5] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[2][0] = 20.;
   MuScleFitUtils::massWindowHalfWidth[2][1] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[2][2] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[2][3] = 0.35;
   MuScleFitUtils::massWindowHalfWidth[2][4] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[2][5] = 0.2;
 
   muonSelector_.reset(new MuScleFitMuonSelector(theMuonLabel_,
                                                 theMuonType_,
                                                 PATmuons_,
                                                 MuScleFitUtils::resfind,
                                                 MuScleFitUtils::speedup,
                                                 genParticlesName_,
                                                 compareToSimTracks_,
                                                 simTracksCollection_,
                                                 MuScleFitUtils::sherpa_,
                                                 debug_));
 
   MuScleFitUtils::backgroundHandler =
       new BackgroundHandler(pset.getParameter<std::vector<int> >("BgrFitType"),
                             pset.getParameter<std::vector<double> >("LeftWindowBorder"),
                             pset.getParameter<std::vector<double> >("RightWindowBorder"),
                             MuScleFitUtils::ResMass,
                             MuScleFitUtils::massWindowHalfWidth[MuScleFitUtils::MuonTypeForCheckMassWindow]);
 
   MuScleFitUtils::crossSectionHandler =
       new CrossSectionHandler(MuScleFitUtils::parCrossSection, MuScleFitUtils::resfind);
 
   // Build cross section scale factors
   // MuScleFitUtils::resfind
 
   MuScleFitUtils::normalizeLikelihoodByEventNumber_ =
       pset.getUntrackedParameter<bool>("NormalizeLikelihoodByEventNumber", true);
   if (debug_ > 0)
     std::cout << "End of MuScleFit constructor" << std::endl;
 
   inputRootTreeFileName_ = pset.getParameter<std::string>("InputRootTreeFileName");
   outputRootTreeFileName_ = pset.getParameter<std::string>("OutputRootTreeFileName");
   maxEventsFromRootTree_ = pset.getParameter<int>("MaxEventsFromRootTree");
 
   MuScleFitUtils::startWithSimplex_ = pset.getParameter<bool>("StartWithSimplex");
   MuScleFitUtils::computeMinosErrors_ = pset.getParameter<bool>("ComputeMinosErrors");
   MuScleFitUtils::minimumShapePlots_ = pset.getParameter<bool>("MinimumShapePlots");
 
   beginOfJobInConstructor();
 }

MuScleFit::~MuScleFit ( )

override

Definition at line 601 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, MuScleFitBase::genMuonPairs_, inputRootTreeFileName_, createfilelist::int, MuScleFitBase::muonPairs_, outputRootTreeFileName_, saveAllToTree_, MuScleFitUtils::SavedPair, findQualityFiles::size, MuScleFitUtils::speedup, MuScleFitBase::theMuonType_, totalEvents_, and RootTreeHandler::writeTree().

                       {
   if (debug_ > 0)
     std::cout << "[MuScleFit]: Destructor" << std::endl;
   std::cout << "Total number of analyzed events = " << totalEvents_ << std::endl;
 
   if (!(outputRootTreeFileName_.empty())) {
     // Save the events to a root tree unless we are reading from the edm root file and the SavedPair size is different from the totalEvents_
     if (!(inputRootTreeFileName_.empty() && (int(MuScleFitUtils::SavedPair.size()) != totalEvents_))) {
       std::cout << "Saving muon pairs to root tree" << std::endl;
       RootTreeHandler rootTreeHandler;
       if (MuScleFitUtils::speedup) {
         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, 0, saveAllToTree_);
         if (debug_ > 0) {
           std::vector<MuonPair>::const_iterator it = muonPairs_.begin();
           std::cout << "[MuScleFit::~MuScleFit] (Destructor)" << std::endl;
           for (; it < muonPairs_.end(); ++it) {
             std::cout << "  Debugging pairs that are going to be written to file" << std::endl;
             std::cout << "  muon1 = " << it->mu1 << std::endl;
             std::cout << "  muon2 = " << it->mu2 << std::endl;
           }
         }
         rootTreeHandler.writeTree(outputRootTreeFileName_, &(muonPairs_), theMuonType_, nullptr, saveAllToTree_);
       } else {
         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, &(MuScleFitUtils::genPair), saveAllToTree_ );
         rootTreeHandler.writeTree(
             outputRootTreeFileName_, &(muonPairs_), theMuonType_, &(genMuonPairs_), saveAllToTree_);
       }
     } else {
       std::cout << "ERROR: events in the vector = " << MuScleFitUtils::SavedPair.size()
                 << " != totalEvents = " << totalEvents_ << std::endl;
     }
   }
 }

Member Function Documentation

void MuScleFit::applyBias	(	reco::Particle::LorentzVector &	mu,
		const int	charge
	)

protected

Apply the bias if needed using the function in MuScleFitUtils.

Definition at line 1396 of file MuScleFit.cc.

References MuScleFitUtils::applyBias(), MuScleFitUtils::BiasType, gather_cfg::cout, MuScleFitBase::debug_, and MuScleFitUtils::goodmuon.

Referenced by selectMuons().

                                                                          {
   if (MuScleFitUtils::BiasType > 0) {
     mu = MuScleFitUtils::applyBias(mu, charge);
     if (debug_ > 0)
       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after bias      Pt = " << mu.Pt() << std::endl;
   }
 }

void MuScleFit::applySmearing ( reco::Particle::LorentzVector & mu )

protected

Apply the smearing if needed using the function in MuScleFitUtils.

Definition at line 1388 of file MuScleFit.cc.

References MuScleFitUtils::applySmearing(), gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::goodmuon, and MuScleFitUtils::SmearType.

Referenced by selectMuons().

                                                            {
   if (MuScleFitUtils::SmearType > 0) {
     mu = MuScleFitUtils::applySmearing(mu);
     if (debug_ > 0)
       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after smearing  Pt = " << mu.Pt() << std::endl;
   }
 }

void MuScleFit::beginOfJobInConstructor ( )

Definition at line 637 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitPlotter::debug, MuScleFitBase::debug_, mps_fire::i, maxLoopNumber, plotter, MuScleFitBase::readProbabilityDistributionsFromFile(), CSCSkim_cfi::rootFileName, contentValuesCheck::ss, AlCaHLTBitMon_QueryRunRegistry::string, MuScleFitBase::theFiles_, MuScleFitBase::theGenInfoRootFileName_, MuScleFitBase::theRootFileName_, and MuScleFitUtils::useProbsFile_.

Referenced by MuScleFit().

 {
   if (debug_ > 0)
     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
   //if(maxLoopNumber>1)
   if (MuScleFitUtils::useProbsFile_) {
     readProbabilityDistributionsFromFile();
   }
 
   if (debug_ > 0)
     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
 
   // Create the root file
   // --------------------
   for (unsigned int i = 0; i < (maxLoopNumber); i++) {
     std::stringstream ss;
     ss << i;
     std::string rootFileName = ss.str() + "_" + theRootFileName_;
     theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE"));
   }
   if (debug_ > 0)
     std::cout << "[MuScleFit]: Root file created" << std::endl;
 
   std::cout << "creating plotter" << std::endl;
   plotter = new MuScleFitPlotter(theGenInfoRootFileName_);
   plotter->debug = debug_;
 }

bool MuScleFit::checkDeltaR	(	reco::Particle::LorentzVector &	genMu,
		reco::Particle::LorentzVector &	recMu
	)

protected

Check if two lorentzVector are near in deltaR.

Definition at line 1355 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::deltaPhi(), PbPb_ZMuSkimMuonDPG_cff::deltaR, and mathSSE::sqrt().

Referenced by duringFastLoop().

                                                                                               {
   //first is always mu-, second is always mu+
   double deltaR =
       sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
            ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
   if (deltaR < 0.01)
     return true;
   else if (debug_ > 0) {
     std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
               << " DOES NOT MATCH with generated muon from resonance: " << std::endl
               << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
   }
   return false;
 }

void MuScleFit::checkParameters ( )

protected

Simple method to check parameters consistency. It aborts the job if the parameters are not consistent.

Definition at line 1406 of file MuScleFit.cc.

Referenced by MuScleFit().

                                 {
   // Fits selection dimension check
   if (MuScleFitUtils::doResolFit.size() != maxLoopNumber) {
     std::cout << "[MuScleFit-Constructor]: wrong size of resolution fits selector = "
               << MuScleFitUtils::doResolFit.size() << std::endl;
     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
     abort();
   }
   if (MuScleFitUtils::doScaleFit.size() != maxLoopNumber) {
     std::cout << "[MuScleFit-Constructor]: wrong size of scale fits selector = " << MuScleFitUtils::doScaleFit.size()
               << std::endl;
     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
     abort();
   }
   if (MuScleFitUtils::doCrossSectionFit.size() != maxLoopNumber) {
     std::cout << "[MuScleFit-Constructor]: wrong size of cross section fits selector = "
               << MuScleFitUtils::doCrossSectionFit.size() << std::endl;
     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
     abort();
   }
   if (MuScleFitUtils::doBackgroundFit.size() != maxLoopNumber) {
     std::cout << "[MuScleFit-Constructor]: wrong size of background fits selector = "
               << MuScleFitUtils::doBackgroundFit.size() << std::endl;
     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
     abort();
   }
 
   // Bias parameters: dimension check
   // --------------------------------
   if ((MuScleFitUtils::BiasType == 1 && MuScleFitUtils::parBias.size() != 2) ||  // linear in pt
       (MuScleFitUtils::BiasType == 2 && MuScleFitUtils::parBias.size() != 2) ||  // linear in |eta|
       (MuScleFitUtils::BiasType == 3 && MuScleFitUtils::parBias.size() != 4) ||  // sinusoidal in phi
       (MuScleFitUtils::BiasType == 4 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and |eta|
       (MuScleFitUtils::BiasType == 5 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and sinusoidal in phi
       (MuScleFitUtils::BiasType == 6 &&
        MuScleFitUtils::parBias.size() != 3) ||  // linear in |eta| and sinusoidal in phi
       (MuScleFitUtils::BiasType == 7 &&
        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and |eta| and sinusoidal in phi
       (MuScleFitUtils::BiasType == 8 && MuScleFitUtils::parBias.size() != 4) ||   // linear in pt and parabolic in |eta|
       (MuScleFitUtils::BiasType == 9 && MuScleFitUtils::parBias.size() != 2) ||   // exponential in pt
       (MuScleFitUtils::BiasType == 10 && MuScleFitUtils::parBias.size() != 3) ||  // parabolic in pt
       (MuScleFitUtils::BiasType == 11 &&
        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and sin in phi with chg
       (MuScleFitUtils::BiasType == 12 &&
        MuScleFitUtils::parBias.size() != 6) ||  // linear in pt and para in plus sin in phi with chg
       (MuScleFitUtils::BiasType == 13 &&
        MuScleFitUtils::parBias.size() != 8) ||  // linear in pt and para in plus sin in phi with chg
       MuScleFitUtils::BiasType < 0 ||
       MuScleFitUtils::BiasType > 13) {
     std::cout << "[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
     abort();
   }
   // Smear parameters: dimension check
   // ---------------------------------
   if ((MuScleFitUtils::SmearType == 1 && MuScleFitUtils::parSmear.size() != 3) ||
       (MuScleFitUtils::SmearType == 2 && MuScleFitUtils::parSmear.size() != 4) ||
       (MuScleFitUtils::SmearType == 3 && MuScleFitUtils::parSmear.size() != 5) ||
       (MuScleFitUtils::SmearType == 4 && MuScleFitUtils::parSmear.size() != 6) ||
       (MuScleFitUtils::SmearType == 5 && MuScleFitUtils::parSmear.size() != 7) ||
       (MuScleFitUtils::SmearType == 6 && MuScleFitUtils::parSmear.size() != 16) ||
       (MuScleFitUtils::SmearType == 7 && !MuScleFitUtils::parSmear.empty()) || MuScleFitUtils::SmearType < 0 ||
       MuScleFitUtils::SmearType > 7) {
     std::cout << "[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
     abort();
   }
   // Protect against bad size of parameters
   // --------------------------------------
   if (MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolFix.size() ||
       MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolOrder.size()) {
     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
     abort();
   }
   if (MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleFix.size() ||
       MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleOrder.size()) {
     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
     abort();
   }
   if (MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionFix.size() ||
       MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionOrder.size()) {
     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
     abort();
   }
   if (MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrFix.size() ||
       MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrOrder.size()) {
     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
     abort();
   }
 
   // Protect against an incorrect number of resonances
   // -------------------------------------------------
   if (MuScleFitUtils::resfind.size() != 6) {
     std::cout << "[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
     abort();
   }
 }

void MuScleFit::duringFastLoop ( )

virtual

This method performs all needed operations on the muon pair. It reads the muons from SavedPair and uses the iev counter to keep track of the event number. The iev is incremented internally and reset to 0 in startingNewLoop.

Definition at line 1080 of file MuScleFit.cc.

References CrossSectionHandler::addParameters(), MuScleFitUtils::applyScale(), checkDeltaR(), compareToSimTracks_, MuScleFitUtils::computeWeight(), gather_cfg::cout, MuScleFitUtils::crossSectionHandler, MuScleFitBase::debug_, MuScleFitUtils::debugMassResol_, MuScleFitUtils::doScaleFit, HLT_2018_cff::eta1, HLT_2018_cff::eta2, fillComparisonHistograms(), dqmdumpme::first, MuScleFitUtils::genPair, mps_fire::i, iev, MuScleFitUtils::iev_, createfilelist::int, dqm-mbProfile::log, loopCounter, MuScleFitBase::mapHisto_, EgHLTOffHistBins_cfi::mass, MuScleFitUtils::massProb(), MuScleFitUtils::massResolComponents, MuScleFitUtils::massResolution(), MuScleFitUtils::parBgr, MuScleFitUtils::parCrossSection, MuScleFitUtils::parResol, MuScleFitUtils::parScale, MuScleFitUtils::parvalue, funct::pow(), TtFullHadEvtBuilder_cfi::prob, HLT_2018_cff::pt1, HLT_2018_cff::pt2, recMu1, recMu2, MuScleFitUtils::ResFound, MuScleFitUtils::resolutionFunctionForVec, MuScleFitUtils::SavedPair, edm::second(), MuScleFitUtils::simPair, findQualityFiles::size, MuScleFitUtils::speedup, and mps_merge::weight.

Referenced by duringLoop(), and endOfLoop().

                                {
   // On loops>0 the two muons are directly obtained from the SavedMuon array
   // -----------------------------------------------------------------------
   MuScleFitUtils::ResFound = false;
   recMu1 = (MuScleFitUtils::SavedPair[iev].first);
   recMu2 = (MuScleFitUtils::SavedPair[iev].second);
 
   //std::cout << "iev = " << iev << ", recMu1 pt = " << recMu1.Pt() << ", recMu2 pt = " << recMu2.Pt() << std::endl;
 
   if (recMu1.Pt() > 0 && recMu2.Pt() > 0) {
     MuScleFitUtils::ResFound = true;
     if (debug_ > 0)
       std::cout << "Ev = " << iev << ": found muons in tree with Pt = " << recMu1.Pt() << " " << recMu2.Pt()
                 << std::endl;
   }
 
   if (debug_ > 0)
     std::cout << "About to start lik par correction and histo filling; ResFound is " << MuScleFitUtils::ResFound
               << std::endl;
   // If resonance found, do the hard work
   // ------------------------------------
   if (MuScleFitUtils::ResFound) {
     // Find weight and reference mass for this muon pair
     // -------------------------------------------------
     // The last parameter = true means that we want to use always the background window to compute the weight,
     // otherwise the probability will be filled only for the resonance region.
     double weight = MuScleFitUtils::computeWeight((recMu1 + recMu2).mass(), iev, true);
     if (debug_ > 0) {
       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": before correction     Pt1 = " << recMu1.Pt()
                 << " Pt2 = " << recMu2.Pt() << std::endl;
     }
     // For successive iterations, correct the muons only if the previous iteration was a scale fit.
     // --------------------------------------------------------------------------------------------
     if (loopCounter > 0) {
       if (MuScleFitUtils::doScaleFit[loopCounter - 1]) {
         recMu1 = (MuScleFitUtils::applyScale(recMu1, MuScleFitUtils::parvalue[loopCounter - 1], -1));
         recMu2 = (MuScleFitUtils::applyScale(recMu2, MuScleFitUtils::parvalue[loopCounter - 1], 1));
       }
     }
     if (debug_ > 0) {
       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": after correction      Pt1 = " << recMu1.Pt()
                 << " Pt2 = " << recMu2.Pt() << std::endl;
     }
 
     reco::Particle::LorentzVector bestRecRes(recMu1 + recMu2);
 
     //Fill histograms
     //------------------
 
     mapHisto_["hRecBestMu"]->Fill(recMu1, -1, weight);
     mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
     mapHisto_["hRecBestMu"]->Fill(recMu2, +1, weight);
     mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
     mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
     // Reconstructed resonance
     mapHisto_["hRecBestRes"]->Fill(bestRecRes, +1, weight);
     mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, +1, 1.);
     //     // Fill histogram of Res mass vs muon variables
     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu1, bestRecRes, -1);
     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu2, bestRecRes, +1);
     //     // Fill also the mass mu+/mu- comparisons
     //     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes);
 
     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, bestRecRes, -1, weight);
     mapHisto_["hRecBestResVSMu"]->Fill(recMu2, bestRecRes, +1, weight);
     // Fill also the mass mu+/mu- comparisons
     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes, weight);
 
     //-- rc 2010 filling histograms for mu+ /mu- ------
     //  mapHisto_["hRecBestResVSMuMinus"]->Fill (recMu1, bestRecRes, -1);
     // mapHisto_["hRecBestResVSMuPlus"]->Fill (recMu2, bestRecRes, +1);
 
     //-- rc 2010 filling histograms MassVsMuEtaPhi------
     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu1, bestRecRes,-1);
     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu2, bestRecRes,+1);
 
     // Fill histogram of Res mass vs Res variables
     // mapHisto_["hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
     mapHisto_["hRecBestResVSRes"]->Fill(bestRecRes, bestRecRes, +1, weight);
 
     std::vector<double>* parval;
     std::vector<double> initpar;
     // Store a pointer to the vector of parameters of the last iteration, or the initial
     // parameters if this is the first iteration
     if (loopCounter == 0) {
       initpar = MuScleFitUtils::parResol;
       initpar.insert(initpar.end(), MuScleFitUtils::parScale.begin(), MuScleFitUtils::parScale.end());
       initpar.insert(initpar.end(), MuScleFitUtils::parCrossSection.begin(), MuScleFitUtils::parCrossSection.end());
       initpar.insert(initpar.end(), MuScleFitUtils::parBgr.begin(), MuScleFitUtils::parBgr.end());
       parval = &initpar;
     } else {
       parval = &(MuScleFitUtils::parvalue[loopCounter - 1]);
     }
 
     //Compute pt resolution w.r.t generated and simulated muons
     //--------------------------------------------------------
     if (!MuScleFitUtils::speedup) {
       //first is always mu-, second is always mu+
       if (checkDeltaR(MuScleFitUtils::genPair[iev].first, recMu1)) {
         fillComparisonHistograms(MuScleFitUtils::genPair[iev].first, recMu1, "Gen", -1);
       }
       if (checkDeltaR(MuScleFitUtils::genPair[iev].second, recMu2)) {
         fillComparisonHistograms(MuScleFitUtils::genPair[iev].second, recMu2, "Gen", +1);
       }
       if (compareToSimTracks_) {
         //first is always mu-, second is always mu+
         if (checkDeltaR(MuScleFitUtils::simPair[iev].first, recMu1)) {
           fillComparisonHistograms(MuScleFitUtils::simPair[iev].first, recMu1, "Sim", -1);
         }
         if (checkDeltaR(MuScleFitUtils::simPair[iev].second, recMu2)) {
           fillComparisonHistograms(MuScleFitUtils::simPair[iev].second, recMu2, "Sim", +1);
         }
       }
     }
 
     // ATTENTION: this was done only when a matching was found. Moved it outside because, genInfo or not, we still want to see the resolution function
     // Fill also the resolution histogramsm using the resolution functions:
     // the parameters are those from the last iteration, as the muons up to this point have also the corrections from the same iteration.
     // Need to use a different array (ForVec), containing functors able to operate on std::vector<double>
     mapHisto_["hFunctionResolPt"]->Fill(
         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), *parval), -1);
     mapHisto_["hFunctionResolCotgTheta"]->Fill(
         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), *parval), -1);
     mapHisto_["hFunctionResolPhi"]->Fill(
         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), *parval), -1);
     mapHisto_["hFunctionResolPt"]->Fill(
         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), *parval), +1);
     mapHisto_["hFunctionResolCotgTheta"]->Fill(
         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), *parval), +1);
     mapHisto_["hFunctionResolPhi"]->Fill(
         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), *parval), +1);
 
     // Compute likelihood histograms
     // -----------------------------
     if (debug_ > 0)
       std::cout << "mass = " << bestRecRes.mass() << std::endl;
     if (weight != 0.) {
       double massResol;
       double prob;
       double deltalike;
       if (loopCounter == 0) {
         std::vector<double> initpar;
         for (int i = 0; i < (int)(MuScleFitUtils::parResol.size()); i++) {
           initpar.push_back(MuScleFitUtils::parResol[i]);
         }
         for (int i = 0; i < (int)(MuScleFitUtils::parScale.size()); i++) {
           initpar.push_back(MuScleFitUtils::parScale[i]);
         }
         //  for (int i=0; i<(int)(MuScleFitUtils::parCrossSection.size()); i++) {
         //    initpar.push_back(MuScleFitUtils::parCrossSection[i]);
         //  }
         MuScleFitUtils::crossSectionHandler->addParameters(initpar);
 
         for (int i = 0; i < (int)(MuScleFitUtils::parBgr.size()); i++) {
           initpar.push_back(MuScleFitUtils::parBgr[i]);
         }
         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, initpar);
         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(), massResol, initpar, true );
         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
                                         bestRecRes.Eta(),
                                         bestRecRes.Rapidity(),
                                         massResol,
                                         initpar,
                                         true,
                                         recMu1.eta(),
                                         recMu2.eta());
       } else {
         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parvalue[loopCounter - 1]);
         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(),
         //                                       massResol, MuScleFitUtils::parvalue[loopCounter-1], true );
         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
                                         bestRecRes.Eta(),
                                         bestRecRes.Rapidity(),
                                         massResol,
                                         MuScleFitUtils::parvalue[loopCounter - 1],
                                         true,
                                         recMu1.eta(),
                                         recMu2.eta());
       }
       if (debug_ > 0)
         std::cout << "inside weight: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
       if (prob > 0) {
         if (debug_ > 0)
           std::cout << "inside prob: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
 
         deltalike = log(prob) * weight;  // NB maximum likelihood --> deltalike is maximized
         mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike);
         mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike);
         mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike);
         mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike);
 
         double recoMass = (recMu1 + recMu2).mass();
         if (recoMass != 0) {
           // IMPORTANT: massResol is not a relative resolution
           mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1);
           mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1);
           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol / recoMass, -1);
           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol / recoMass, +1);
         }
 
         if (MuScleFitUtils::debugMassResol_) {
           mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1);
           mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1);
           mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1);
           mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1);
           mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1);
           mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1);
         }
 
         if (!MuScleFitUtils::speedup) {
           double genMass = (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second).mass();
           // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
           if (genMass != 0) {
             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].first),
                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
                                            -1);
             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].second),
                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
                                            +1);
             double diffMass = (recoMass - genMass) / genMass;
             // double diffMass = recoMass - genMass;
             // Fill if for both muons
             double pt1 = recMu1.pt();
             double eta1 = recMu1.eta();
             double pt2 = recMu2.pt();
             double eta2 = recMu2.eta();
             // This is to avoid nan
             if (diffMass == diffMass) {
               // Mass relative difference vs Pt and Eta. To be used to extract the true mass resolution
               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
               // This is used for the covariance comparison
               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
             } else {
               std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
             }
           }
           // Fill with mass resolution from resolution function
           double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
           mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
           mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
         }
 
         mapHisto_["hMass_P"]->Fill(bestRecRes.mass(), prob);
         if (debug_ > 0)
           std::cout << "mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
         mapHisto_["hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
 
         mapHisto_["hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
         mapHisto_["hMassProbVsMu"]->Fill(recMu1, bestRecRes, -1, prob);
         mapHisto_["hMassProbVsMu"]->Fill(recMu2, bestRecRes, +1, prob);
         mapHisto_["hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu1, bestRecRes, -1, prob);
         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu2, bestRecRes, +1, prob);
       }
     }
   }  // end if ResFound
 
   // Fill the pair
   // -------------
   if (loopCounter > 0) {
     if (debug_ > 0)
       std::cout << "[MuScleFit]: filling the pair" << std::endl;
     MuScleFitUtils::SavedPair[iev] = std::make_pair(recMu1, recMu2);
   }
 
   iev++;
   MuScleFitUtils::iev_++;
 
   // return kContinue;
 }

edm::EDLooper::Status MuScleFit::duringLoop	(	const edm::Event &	,
		const edm::EventSetup &
	)

overridevirtual

Called after all event modules have had a chance to process the edm::Event.

Implements edm::EDLooper.

Definition at line 782 of file MuScleFit.cc.

References edm::HLTGlobalStatus::accept(), gather_cfg::cout, MuScleFitBase::debug_, duringFastLoop(), relativeConstraints::empty, fastLoop, edm::Event::getRun(), HltComparatorCreateWorkflow::hltConfig, mps_fire::i, HLTConfigProvider::init(), inputRootTreeFileName_, edm::EDLooperBase::kContinue, loopCounter, negateTrigger_, selectMuons(), edm::TriggerNames::size(), findQualityFiles::size, AlCaHLTBitMon_QueryRunRegistry::string, totalEvents_, HLTConfigProvider::triggerIndex(), edm::TriggerNames::triggerName(), L1TEGammaOffline_cfi::triggerNames, triggerPath_, triggerResults, triggerResultsLabel_, and triggerResultsProcess_.

                                                                                               {
   edm::Handle<edm::TriggerResults> triggerResults;
   event.getByLabel(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()), triggerResults);
   //event.getByLabel(InputTag(triggerResultsLabel_),triggerResults);
   bool isFired = false;
 
   if (triggerPath_[0].empty())
     isFired = true;
   else if (triggerPath_[0] == "All") {
     isFired = triggerResults->accept();
     if (debug_ > 0)
       std::cout << "Trigger " << isFired << std::endl;
   } else {
     bool changed;
     HLTConfigProvider hltConfig;
     hltConfig.init(event.getRun(), eventSetup, triggerResultsProcess_, changed);
 
     const edm::TriggerNames& triggerNames = event.triggerNames(*triggerResults);
 
     for (unsigned i = 0; i < triggerNames.size(); i++) {
       const std::string& hltName = triggerNames.triggerName(i);
 
       // match the path in the pset with the true name of the trigger
       for (unsigned int ipath = 0; ipath < triggerPath_.size(); ipath++) {
         if (hltName.find(triggerPath_[ipath]) != std::string::npos) {
           unsigned int triggerIndex(hltConfig.triggerIndex(hltName));
 
           // triggerIndex must be less than the size of HLTR or you get a CMSException: _M_range_check
           if (triggerIndex < triggerResults->size()) {
             isFired = triggerResults->accept(triggerIndex);
             if (debug_ > 0)
               std::cout << triggerPath_[ipath] << " " << hltName << " " << isFired << std::endl;
           }
         }  // end if (matching the path in the pset with the true trigger name
       }
     }
   }
 
   if (negateTrigger_ && isFired)
     return kContinue;
   else if (!(negateTrigger_) && !isFired)
     return kContinue;
 
 #ifdef USE_CALLGRIND
   CALLGRIND_START_INSTRUMENTATION;
 #endif
 
   if (debug_ > 0) {
     std::cout << "[MuScleFit-duringLoop]: loopCounter = " << loopCounter << " Run: " << event.id().run()
               << " Event: " << event.id().event() << std::endl;
   }
 
   // On the first iteration we read the bank, otherwise we fetch the information from the muon tree
   // ------------------------------------ Important Note --------------------------------------- //
   // The fillMuonCollection method applies any smearing or bias to the muons, so we NEVER use
   // unbiased muons.
   // ----------------------------------------------------------------------------------------------
   if (loopCounter == 0) {
     if (!fastLoop || inputRootTreeFileName_.empty()) {
       if (debug_ > 0)
         std::cout << "Reading from edm event" << std::endl;
       selectMuons(event);
       duringFastLoop();
       ++totalEvents_;
     }
   }
 
   return kContinue;
 
 #ifdef USE_CALLGRIND
   CALLGRIND_STOP_INSTRUMENTATION;
   CALLGRIND_DUMP_STATS;
 #endif
 }

void MuScleFit::endOfFastLoop ( const unsigned int iLoop )

virtual

Definition at line 747 of file MuScleFit.cc.

References MuScleFitBase::clearHistoMap(), gather_cfg::cout, loopCounter, MuScleFitUtils::minimizeLikelihood(), plotter, MuScleFitBase::theFiles_, and MuScleFitBase::writeHistoMap().

Referenced by endOfLoop().

                                                       {
   // std::cout<< "Inside endOfFastLoop, iLoop = " << iLoop << " and loopCounter = " << loopCounter << std::endl;
 
   if (loopCounter == 0) {
     // plotter->writeHistoMap();
     // The destructor will call the writeHistoMap after the cd to the output file
     delete plotter;
   }
 
   std::cout << "Ending loop # " << iLoop << std::endl;
 
   // Write the histos to file
   // ------------------------
   // theFiles_[iLoop]->cd();
   writeHistoMap(iLoop);
 
   // Likelihood minimization to compute corrections
   // ----------------------------------------------
   // theFiles_[iLoop]->cd();
   TDirectory* likelihoodDir = theFiles_[iLoop]->mkdir("likelihood");
   likelihoodDir->cd();
   MuScleFitUtils::minimizeLikelihood();
 
   // ATTENTION, this was put BEFORE the minimizeLikelihood. Check for problems.
   theFiles_[iLoop]->Close();
   // ATTENTION: Check that this delete does not give any problem
   delete theFiles_[iLoop];
 
   // Clear the histos
   // ----------------
   clearHistoMap();
 }

void MuScleFit::endOfJob ( )

overridevirtual

Reimplemented from edm::EDLooperBase.

Definition at line 669 of file MuScleFit.cc.

References gather_cfg::cout, and MuScleFitBase::debug_.

                          {
   if (debug_ > 0)
     std::cout << "[MuScleFit]: endOfJob" << std::endl;
 }

edm::EDLooper::Status MuScleFit::endOfLoop	(	const edm::EventSetup &	eventSetup,
		unsigned int	iLoop
	)

override

Definition at line 703 of file MuScleFit.cc.

References gather_cfg::cout, duringFastLoop(), endOfFastLoop(), fastLoop, iev, inputRootTreeFileName_, edm::EDLooperBase::kContinue, edm::EDLooperBase::kStop, maxEventsFromRootTree_, maxLoopNumber, MuScleFitUtils::SavedPair, selectMuons(), startingNewLoop(), and totalEvents_.

                                                                                           {
   unsigned int iFastLoop = 1;
 
   // Read the events from the root tree if requested
   if (!(inputRootTreeFileName_.empty())) {
     selectMuons(maxEventsFromRootTree_, inputRootTreeFileName_);
     // When reading from local file all the loops are done here
     totalEvents_ = MuScleFitUtils::SavedPair.size();
     iFastLoop = 0;
   } else {
     endOfFastLoop(iLoop);
   }
 
   // If a fastLoop is required we do all the remaining iterations here
   if (fastLoop == true) {
     for (; iFastLoop < maxLoopNumber; ++iFastLoop) {
       std::cout << "Starting fast loop number " << iFastLoop << std::endl;
 
       // In the first loop is called by the framework
       // if( iFastLoop > 0 ) {
       startingNewLoop(iFastLoop);
       // }
 
       // std::vector<std::pair<lorentzVector,lorentzVector> >::const_iterator it = MuScleFitUtils::SavedPair.begin();
       // for( ; it != SavedPair.end(); ++it ) {
       while (iev < totalEvents_) {
         if (iev % 50000 == 0) {
           std::cout << "Fast looping on event number " << iev << std::endl;
         }
         // This reads muons from SavedPair using iev to keep track of the event
         duringFastLoop();
       }
       std::cout << "End of fast loop number " << iFastLoop << ". Ran on " << iev << " events" << std::endl;
       endOfFastLoop(iFastLoop);
     }
   }
 
   if (iFastLoop >= maxLoopNumber - 1) {
     return kStop;
   } else {
     return kContinue;
   }
 }

void MuScleFit::fillComparisonHistograms	(	const reco::Particle::LorentzVector &	genMu,
		const reco::Particle::LorentzVector &	recoMu,
		const std::string &	inputName,
		const int	charge
	)

protected

Fill the reco vs gen and reco vs sim comparison histograms.

Definition at line 1370 of file MuScleFit.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, funct::cos(), MuScleFitUtils::deltaPhiNoFabs(), MuScleFitBase::mapHisto_, Skims_PA_cff::name, funct::sin(), and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by duringFastLoop().

                                                            {
   std::string name(inputName + "VSMu");
   mapHisto_["hResolPt" + name]->Fill(recMu, (-genMu.Pt() + recMu.Pt()) / genMu.Pt(), charge);
   mapHisto_["hResolTheta" + name]->Fill(recMu, (-genMu.Theta() + recMu.Theta()), charge);
   mapHisto_["hResolCotgTheta" + name]->Fill(
       recMu, (-cos(genMu.Theta()) / sin(genMu.Theta()) + cos(recMu.Theta()) / sin(recMu.Theta())), charge);
   mapHisto_["hResolEta" + name]->Fill(recMu, (-genMu.Eta() + recMu.Eta()), charge);
   mapHisto_["hResolPhi" + name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);
 
   // Fill only if it was matched to a genMu and this muon is valid
   if ((genMu.Pt() != 0) && (recMu.Pt() != 0)) {
     mapHisto_["hPtRecoVsPt" + inputName]->Fill(genMu.Pt(), recMu.Pt());
   }
 }

template<typename T >

std::vector< MuScleFitMuon > MuScleFit::fillMuonCollection ( const std::vector< T > & tracks )

Definition at line 314 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitUtils::goodmuon, MuScleFitUtils::mMu2, amptDefaultParameters_cff::mu, HLT_2018_cff::muon, PDWG_BPHSkim_cff::muons, mathSSE::sqrt(), and HLT_2018_cff::track.

                                                                                  {
   std::vector<MuScleFitMuon> muons;
   typename std::vector<T>::const_iterator track;
   for (track = tracks.begin(); track != tracks.end(); ++track) {
     reco::Particle::LorentzVector mu;
     mu = reco::Particle::LorentzVector(
         track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
     // Apply smearing if needed, and then bias
     // ---------------------------------------
     MuScleFitUtils::goodmuon++;
     if (debug_ > 0)
       std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon << ": initial value   Pt = " << mu.Pt()
                 << std::endl;
 
     applySmearing(mu);
     applyBias(mu, track->charge());
     if (debug_ > 0)
       std::cout << "track charge: " << track->charge() << std::endl;
 
     Double_t hitsTk = track->innerTrack()->hitPattern().numberOfValidTrackerHits();
     Double_t hitsMuon = track->innerTrack()->hitPattern().numberOfValidMuonHits();
     Double_t ptError = track->innerTrack()->ptError();
     MuScleFitMuon muon(mu, track->charge(), ptError, hitsTk, hitsMuon, false);
     if (debug_ > 0) {
       std::cout << "[MuScleFit::fillMuonCollection]" << std::endl;
       std::cout << "  muon = " << muon << std::endl;
     }
 
     // Store modified muon
     // -------------------
     muons.push_back(muon);
   }
   return muons;
 }

void MuScleFit::selectMuons ( const edm::Event & event )

protected

Selects the muon pairs and fills the SavedPair and (if needed) the genPair vector. This version reads the events from the edm root file and performs a selection of the muons according to the parameters in the cfg.

Definition at line 857 of file MuScleFit.cc.

References L1TStage2uGTEmulatorClient_cff::BX, gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::findBestRecoRes(), dqmdumpme::first, MuScleFitBase::genMuonPairs_, MuScleFitUtils::genPair, edm::HandleBase::isValid(), MuScleFitBase::muonPairs_, PDWG_BPHSkim_cff::muons, muonSelector_, plotter, puInfoSrc_, recMu1, recMu2, recMuScleMu1, recMuScleMu2, MuScleFitUtils::ResFound, edm::Event::run(), MuScleFitUtils::SavedPair, MuScleFitUtils::SavedPairMuScleFitMuons, edm::second(), MuScleFitUtils::simPair, MuScleFitUtils::speedup, vertexSrc_, pwdgSkimBPark_cfi::vertices, and GenEventInfoProduct::weight().

Referenced by duringLoop(), and endOfLoop().

                                                  {
   recMu1 = reco::Particle::LorentzVector(0, 0, 0, 0);
   recMu2 = reco::Particle::LorentzVector(0, 0, 0, 0);
 
   std::vector<MuScleFitMuon> muons;
   muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
   //  plotter->fillRec(muons); // @EM method already invoked inside MuScleFitMuonSelector::selectMuons()
 
   if (debug_ > 0) {
     std::cout << "[MuScleFit::selectMuons] Debugging muons collections after call to muonSelector_->selectMuons"
               << std::endl;
     int iMu = 0;
     for (std::vector<MuScleFitMuon>::const_iterator it = muons.begin(); it < muons.end(); ++it) {
       std::cout << "  - muon n. " << iMu << " = " << (*it) << std::endl;
       ++iMu;
     }
   }
 
   // Find the two muons from the resonance, and set ResFound bool
   // ------------------------------------------------------------
   std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes = MuScleFitUtils::findBestRecoRes(muons);
 
   if (MuScleFitUtils::ResFound) {
     if (debug_ > 0) {
       std::cout << std::setprecision(9) << "Pt after findbestrecores: " << (recMuFromBestRes.first).Pt() << " "
                 << (recMuFromBestRes.second).Pt() << std::endl;
       std::cout << "recMu1 = " << recMu1 << std::endl;
       std::cout << "recMu2 = " << recMu2 << std::endl;
     }
     recMu1 = recMuFromBestRes.first.p4();
     recMu2 = recMuFromBestRes.second.p4();
     recMuScleMu1 = recMuFromBestRes.first;
     recMuScleMu2 = recMuFromBestRes.second;
 
     if (debug_ > 0) {
       std::cout << "after recMu1 = " << recMu1 << std::endl;
       std::cout << "after recMu2 = " << recMu2 << std::endl;
       std::cout << "mu1.pt = " << recMu1.Pt() << std::endl;
       std::cout << "mu2.pt = " << recMu2.Pt() << std::endl;
       std::cout << "after recMuScleMu1 = " << recMuScleMu1 << std::endl;
       std::cout << "after recMuScleMu2 = " << recMuScleMu2 << std::endl;
     }
     MuScleFitUtils::SavedPair.push_back(std::make_pair(recMu1, recMu2));
     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(recMuScleMu1, recMuScleMu2));
   } else {
     MuScleFitUtils::SavedPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(MuScleFitMuon(), MuScleFitMuon()));
   }
   // Save the events also in the external tree so that it can be saved late
 
   // Fetch extra information (per event)
   UInt_t the_NVtx(0);
   Int_t the_numPUvtx(0);
   Float_t the_TrueNumInteractions(0);
 
   // Fill pile-up related informations
   // --------------------------------
   edm::Handle<std::vector<PileupSummaryInfo> > puInfo;
   event.getByLabel(puInfoSrc_, puInfo);
   if (puInfo.isValid()) {
     std::vector<PileupSummaryInfo>::const_iterator PVI;
     for (PVI = puInfo->begin(); PVI != puInfo->end(); ++PVI) {
       int BX = PVI->getBunchCrossing();
       if (BX == 0) {  // "0" is the in-time crossing, negative values are the early crossings, positive are late
         the_TrueNumInteractions = PVI->getTrueNumInteractions();
         the_numPUvtx = PVI->getPU_NumInteractions();
       }
     }
   }
 
   edm::Handle<std::vector<reco::Vertex> > vertices;
   event.getByLabel(vertexSrc_, vertices);
   if (vertices.isValid()) {
     std::vector<reco::Vertex>::const_iterator itv;
     // now, count vertices
     for (itv = vertices->begin(); itv != vertices->end(); ++itv) {
       // require that the vertex meets certain criteria
       if (itv->ndof() < 5)
         continue;
       if (fabs(itv->z()) > 50.0)
         continue;
       if (fabs(itv->position().rho()) > 2.0)
         continue;
       ++the_NVtx;
     }
   }
 
   // get the MC event weight
   edm::Handle<GenEventInfoProduct> genEvtInfo;
   event.getByLabel("generator", genEvtInfo);
   double the_genEvtweight = 1.;
   if (genEvtInfo.isValid()) {
     the_genEvtweight = genEvtInfo->weight();
   }
 
   muonPairs_.push_back(MuonPair(
       MuScleFitUtils::SavedPairMuScleFitMuons.back().first,
       MuScleFitUtils::SavedPairMuScleFitMuons.back().second,
       MuScleFitEvent(
           event.run(), event.id().event(), the_genEvtweight, the_numPUvtx, the_TrueNumInteractions, the_NVtx)));
   // Fill the internal genPair tree from the external one
   if (MuScleFitUtils::speedup == false) {
     MuScleFitUtils::genPair.push_back(std::make_pair(genMuonPairs_.back().mu1.p4(), genMuonPairs_.back().mu2.p4()));
   }
 }

void MuScleFit::selectMuons	(	const int	maxEvents,
		const TString &	treeFileName
	)

protected

Selects the muon pairs and fills the SavedPair and (if needed) the genPair vector. This version reads the events from a tree in the file specified in the cfg. The tree only contains one muon pair per event. This means that no selection is performed and we use preselected muons.

Definition at line 963 of file MuScleFit.cc.

References applyBias(), applySmearing(), MuScleFitUtils::BiasType, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, MuScleFitUtils::deltaPhi(), MuScleFitUtils::deltaPhiMaxCut_, MuScleFitUtils::deltaPhiMinCut_, HLT_2018_cff::eta1, HLT_2018_cff::eta2, MuScleFitPlotter::fillTreeGen(), MuScleFitPlotter::fillTreeRec(), MuScleFitBase::genMuonPairs_, MuScleFitUtils::genMuscleFitPair, MuScleFitUtils::genPair, inputRootTreeFileName_, MuScleFitUtils::maxMuonEtaFirstRange_, MuScleFitUtils::maxMuonEtaSecondRange_, MuScleFitUtils::maxMuonPt_, MuScleFitUtils::minMuonEtaFirstRange_, MuScleFitUtils::minMuonEtaSecondRange_, MuScleFitUtils::minMuonPt_, MuScleFitBase::muonPairs_, plotter, HLT_2018_cff::pt1, HLT_2018_cff::pt2, RootTreeHandler::readTree(), MuScleFitUtils::SavedPair, MuScleFitUtils::SavedPairMuScleFitMuons, MuScleFitUtils::separateRanges_, MuScleFitUtils::SmearType, MuScleFitUtils::speedup, and MuScleFitBase::theMuonType_.

                                                                             {
   std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
   RootTreeHandler rootTreeHandler;
   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
   if (MuScleFitUtils::speedup) {
     rootTreeHandler.readTree(
         maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPairMuScleFitMuons), theMuonType_, &evtRun);
   } else {
     rootTreeHandler.readTree(maxEvents,
                              inputRootTreeFileName_,
                              &(MuScleFitUtils::SavedPairMuScleFitMuons),
                              theMuonType_,
                              &evtRun,
                              &(MuScleFitUtils::genMuscleFitPair));
   }
   // Now loop on all the pairs and apply any smearing and bias if needed
   std::vector<std::pair<unsigned int, unsigned long long> >::iterator evtRunIt = evtRun.begin();
   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator it = MuScleFitUtils::SavedPairMuScleFitMuons.begin();
   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator genIt;
   if (MuScleFitUtils::speedup == false)
     genIt = MuScleFitUtils::genMuscleFitPair.begin();
   for (; it != MuScleFitUtils::SavedPairMuScleFitMuons.end(); ++it, ++evtRunIt) {
     // Apply any cut if requested
     // Note that cuts here are only applied to already selected muons. They should not be used unless
     // you are sure that the difference is negligible (e.g. the number of events with > 2 muons is negligible).
     double pt1 = it->first.pt();
     //std::cout << "pt1 = " << pt1 << std::endl;
     double pt2 = it->second.pt();
     //std::cout << "pt2 = " << pt2 << std::endl;
     double eta1 = it->first.eta();
     //std::cout << "eta1 = " << eta1 << std::endl;
     double eta2 = it->second.eta();
     //std::cout << "eta2 = " << eta2 << std::endl;
     // If they don't pass the cuts set to null vectors
     bool dontPass = false;
     bool eta1InFirstRange;
     bool eta2InFirstRange;
     bool eta1InSecondRange;
     bool eta2InSecondRange;
 
     int ch1 = it->first.charge();
     int ch2 = it->second.charge();
 
     if (MuScleFitUtils::separateRanges_) {
       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
       eta1InSecondRange =
           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
       eta2InSecondRange =
           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
 
       // This is my logic, which should be erroneous, but certainly simpler...
       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
             ((eta1InFirstRange && eta2InSecondRange && ch1 >= ch2) ||
              (eta1InSecondRange && eta2InFirstRange && ch1 < ch2))))
         dontPass = true;
     } else {
       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
       eta1InSecondRange =
           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
       eta2InSecondRange =
           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
             (((eta1InFirstRange && !eta2InFirstRange) && (eta2InSecondRange && !eta1InSecondRange) && ch1 >= ch2) ||
              ((eta2InFirstRange && !eta1InFirstRange) && (eta1InSecondRange && !eta2InSecondRange) && ch1 < ch2))))
         dontPass = true;
     }
 
     // Additional check on deltaPhi
     double deltaPhi = MuScleFitUtils::deltaPhi(it->first.phi(), it->second.phi());
     if ((deltaPhi <= MuScleFitUtils::deltaPhiMinCut_) || (deltaPhi >= MuScleFitUtils::deltaPhiMaxCut_))
       dontPass = true;
 
     lorentzVector vec1 = it->first.p4();
     lorentzVector vec2 = it->second.p4();
     if (ch1 >= ch2) {
       lorentzVector vectemp = vec1;
       vec1 = vec2;
       vec2 = vectemp;
     }
 
     if (!dontPass) {
       // First is always mu-, second mu+
       if ((MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0)) {
         applySmearing(vec1);
         applyBias(vec1, -1);
         applySmearing(vec2);
         applyBias(vec2, 1);
       }
 
       MuScleFitUtils::SavedPair.push_back(std::make_pair(vec1, vec2));
     }
 
     //FIXME: we loose the additional information besides the 4-momenta
     muonPairs_.push_back(MuonPair(
         MuScleFitMuon(vec1, -1),
         MuScleFitMuon(vec2, +1),
         MuScleFitEvent(
             (*evtRunIt).first, (*evtRunIt).second, 0, 0, 0, 0))  // FIXME: order of event and run number mixed up!
     );
 
     // Fill the internal genPair tree from the external one
     if (!MuScleFitUtils::speedup) {
       MuScleFitUtils::genPair.push_back(std::make_pair(genIt->first.p4(), genIt->second.p4()));
       genMuonPairs_.push_back(GenMuonPair(genIt->first.p4(), genIt->second.p4(), 0));
       ++genIt;
     }
   }
   plotter->fillTreeRec(MuScleFitUtils::SavedPair);
   if (!(MuScleFitUtils::speedup)) {
     plotter->fillTreeGen(MuScleFitUtils::genPair);
   }
 }

bool MuScleFit::selGlobalMuon ( const pat::Muon * aMuon )

protected

Function for onia selections.

Definition at line 1502 of file MuScleFit.cc.

References pat::Muon::globalTrack(), pat::Muon::innerTrack(), pat::Muon::muonID(), reco::HitPattern::numberOfValidMuonHits(), AlCaHLTBitMon_ParallelJobs::p, reco::HitPattern::pixelLayersWithMeasurement(), and data-class-funcs::q.

                                                   {
   reco::TrackRef iTrack = aMuon->innerTrack();
   const reco::HitPattern& p = iTrack->hitPattern();
 
   reco::TrackRef gTrack = aMuon->globalTrack();
   const reco::HitPattern& q = gTrack->hitPattern();
 
   return (  //isMuonInAccept(aMuon) &&// no acceptance cuts!
       iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 && q.numberOfValidMuonHits() > 0 &&
       iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
       fabs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
       fabs(iTrack->dz()) < 15.0);   //should be done w.r.t. PV!
 }

bool MuScleFit::selTrackerMuon ( const pat::Muon * aMuon )

protected

Definition at line 1517 of file MuScleFit.cc.

References DEFINE_FWK_LOOPER, pat::Muon::innerTrack(), pat::Muon::muonID(), AlCaHLTBitMon_ParallelJobs::p, and reco::HitPattern::pixelLayersWithMeasurement().

                                                    {
   reco::TrackRef iTrack = aMuon->innerTrack();
   const reco::HitPattern& p = iTrack->hitPattern();
 
   return (  //isMuonInAccept(aMuon) // no acceptance cuts!
       iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
       fabs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
       fabs(iTrack->dz()) < 15.0);   //should be done w.r.t. PV!
 }

void MuScleFit::startingNewLoop ( unsigned int int )

overridevirtual

Called before system starts to loop over the events. The argument is a count of how many loops have been processed. For the first time through the events the argument will be 0.

Implements edm::EDLooperBase.

Definition at line 676 of file MuScleFit.cc.

References MuScleFitUtils::counter_resprob, gather_cfg::cout, MuScleFitBase::debug_, MuScleFitBase::fillHistoMap(), MuScleFitUtils::goodmuon, iev, MuScleFitUtils::iev_, MuScleFitUtils::loopCounter, loopCounter, MuScleFitUtils::oldNormalization_, and MuScleFitBase::theFiles_.

Referenced by endOfLoop().

                                                   {
   if (debug_ > 0)
     std::cout << "[MuScleFit]: Starting loop # " << iLoop << std::endl;
 
   // Number of muons used
   // --------------------
   MuScleFitUtils::goodmuon = 0;
 
   // Counters for problem std::cout-ing
   // -----------------------------
   MuScleFitUtils::counter_resprob = 0;
 
   // Create the root file
   // --------------------
   fillHistoMap(theFiles_[iLoop], iLoop);
 
   loopCounter = iLoop;
   MuScleFitUtils::loopCounter = loopCounter;
 
   iev = 0;
   MuScleFitUtils::iev_ = 0;
 
   MuScleFitUtils::oldNormalization_ = 0;
 }

template<typename T >

void MuScleFit::takeSelectedMuonType	(	const T &	muon,
		std::vector< reco::Track > &	tracks
	)

protected

Template method used to fill the track collection starting from reco::muons or pat::muons.

Definition at line 350 of file MuScleFit.cc.

                                                                                 {
   // std::cout<<"muon "<<muon->isGlobalMuon()<<muon->isStandAloneMuon()<<muon->isTrackerMuon()<<std::endl;
   //NNBB: one muon can be of many kinds at once but with the theMuonType_ we are sure
   // to avoid double counting of the same muon
   if (muon->isGlobalMuon() && theMuonType_ == 1)
     tracks.push_back(*(muon->globalTrack()));
   else if (muon->isStandAloneMuon() && theMuonType_ == 2)
     tracks.push_back(*(muon->outerTrack()));
   else if (muon->isTrackerMuon() && theMuonType_ == 3)
     tracks.push_back(*(muon->innerTrack()));
 
   else if (theMuonType_ == 10 && !(muon->isStandAloneMuon()))  //particular case!!
     tracks.push_back(*(muon->innerTrack()));
   else if (theMuonType_ == 11 && muon->isGlobalMuon())
     tracks.push_back(*(muon->innerTrack()));
   else if (theMuonType_ == 13 && muon->isTrackerMuon())
     tracks.push_back(*(muon->innerTrack()));
 }