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

void	esPrefetchAsync (WaitingTaskHolder iTask, EventSetupImpl const &iImpl, Transition iTrans, ServiceToken const &iToken) const

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

EDLooperBase &	operator= (EDLooperBase const &)=delete

void	prefetchAsync (WaitingTaskHolder iTask, ServiceToken const &token, Transition iTrans, Principal const &iPrincipal, EventSetupImpl const &iImpl) const

void	prepareForNextLoop (eventsetup::EventSetupProvider *esp)

void	setActionTable (ExceptionToActionTable const *actionTable)

void	setModuleChanger (ModuleChanger *)

	~EDLooperBase () noexcept(false) override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESResolverIndex const *	esGetTokenIndices (edm::Transition iTrans) const

std::vector< ESResolverIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const

std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::array< std::vector< ModuleDescription const > , NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, ProductRegistry const &preg, std::map< std::string, ModuleDescription const *> const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProductResolverIndices const &)

virtual	~EDConsumerBase () 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 Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()

template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)

template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

void	resetItemsToGetFrom (BranchType iType)

Protected Attributes
bool	compareToSimTracks_

bool	fastLoop

edm::EDGetTokenT< GenEventInfoProduct >	genEvtInfoToken_

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_

edm::EDGetTokenT< std::vector< PileupSummaryInfo > >	puInfoToken_

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::EDGetTokenT< edm::TriggerResults >	triggerResultsToken_

edm::InputTag	vertexSrc_

edm::EDGetTokenT< reco::VertexCollection >	vertexToken_

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

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

virtual	~MuScleFitBase () noexcept(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_

int	theCompressionSettings_

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 185 of file MuScleFit.cc.

Constructor & Destructor Documentation

◆ MuScleFit()

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

Definition at line 378 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");
 
   triggerResultsToken_ =
       consumes<edm::TriggerResults>(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()));
 
   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");
   puInfoToken_ = consumes<std::vector<PileupSummaryInfo> >(puInfoSrc_);
 
   vertexSrc_ = pset.getUntrackedParameter<edm::InputTag>("PrimaryVertexCollection");
   vertexToken_ = consumes<reco::VertexCollection>(vertexSrc_);
 
   genEvtInfoToken_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
 
   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;
 
   edm::ConsumesCollector iC = consumesCollector();
   muonSelector_ = std::make_unique<MuScleFitMuonSelector>(iC,
                                                           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::~MuScleFit ( )

override

Definition at line 619 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

◆ applyBias()

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

protected

Apply the bias if needed using the function in MuScleFitUtils.

Definition at line 1414 of file MuScleFit.cc.

References MuScleFitUtils::applyBias(), MuScleFitUtils::BiasType, ALCARECOTkAlJpsiMuMu_cff::charge, gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::goodmuon, and amptDefaultParameters_cff::mu.

Referenced by fillMuonCollection(), and 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;
   }
 }

◆ applySmearing()

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

protected

Apply the smearing if needed using the function in MuScleFitUtils.

Definition at line 1406 of file MuScleFit.cc.

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

Referenced by fillMuonCollection(), and selectMuons().

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

◆ beginOfJobInConstructor()

void MuScleFit::beginOfJobInConstructor ( )

Definition at line 655 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::theCompressionSettings_, 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_;
     if (theCompressionSettings_ > -1) {
       theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE", "", theCompressionSettings_));
     } else {
       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_;
 }

◆ checkDeltaR()

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

protected

Check if two lorentzVector are near in deltaR.

Definition at line 1373 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;
 }

◆ checkParameters()

void MuScleFit::checkParameters ( )

protected

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

Definition at line 1424 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();
   }
 }

◆ duringFastLoop()

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 1097 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_2023v12_cff::eta1, HLT_2023v12_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, conifer::pow(), TtFullHadEvtBuilder_cfi::prob, HLT_2023v12_cff::pt1, HLT_2023v12_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;
         initpar.reserve((int)(MuScleFitUtils::parResol.size()));
         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;
 }

◆ duringLoop()

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 804 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, duringFastLoop(), relativeConstraints::empty, options_cfi::eventSetup, fastLoop, HltComparatorCreateWorkflow::hltConfig, mps_fire::i, inputRootTreeFileName_, edm::EDLooperBase::kContinue, loopCounter, negateTrigger_, selectMuons(), findQualityFiles::size, AlCaHLTBitMon_QueryRunRegistry::string, totalEvents_, L1TEGammaOffline_cfi::triggerNames, triggerPath_, triggerResults, triggerResultsProcess_, and triggerResultsToken_.

                                                                                               {
   edm::Handle<edm::TriggerResults> triggerResults = event.getHandle(triggerResultsToken_);
   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
 }

◆ endOfFastLoop()

void MuScleFit::endOfFastLoop ( const unsigned int iLoop )

virtual

Definition at line 769 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();
 }

◆ endOfJob()

void MuScleFit::endOfJob ( )

overridevirtual

Reimplemented from edm::EDLooperBase.

Definition at line 691 of file MuScleFit.cc.

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

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

◆ endOfLoop()

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

overridevirtual

Called after the system has finished one loop over the events. Thar argument is a count of how many loops have been processed before this loo. For the first time through the events the argument will be 0.

Implements edm::EDLooperBase.

Definition at line 725 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;
   }
 }

◆ fillComparisonHistograms()

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 1388 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());
   }
 }

◆ fillMuonCollection()

template<typename T >

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

Definition at line 321 of file MuScleFit.cc.

References applyBias(), applySmearing(), gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::goodmuon, MuScleFitUtils::mMu2, amptDefaultParameters_cff::mu, HLT_2023v12_cff::muon, PDWG_BPHSkim_cff::muons, mathSSE::sqrt(), HLT_2023v12_cff::track, and pwdgSkimBPark_cfi::tracks.

                                                                                  {
   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;
 }

◆ selectMuons() [1/2]

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 877 of file MuScleFit.cc.

References L1TStage2uGTEmulatorClient_cff::BX, gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::findBestRecoRes(), dqmdumpme::first, genEvtInfoToken_, MuScleFitBase::genMuonPairs_, MuScleFitUtils::genPair, edm::HandleBase::isValid(), MuScleFitBase::muonPairs_, PDWG_BPHSkim_cff::muons, muonSelector_, plotter, puInfoToken_, recMu1, recMu2, recMuScleMu1, recMuScleMu2, MuScleFitUtils::ResFound, MuScleFitUtils::SavedPair, MuScleFitUtils::SavedPairMuScleFitMuons, edm::second(), MuScleFitUtils::simPair, MuScleFitUtils::speedup, vertexToken_, AlignmentTracksFromVertexSelector_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.getHandle(puInfoToken_);
   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.getHandle(vertexToken_);
   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.getHandle(genEvtInfoToken_);
   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()));
   }
 }

◆ selectMuons() [2/2]

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 980 of file MuScleFit.cc.

References applyBias(), applySmearing(), MuScleFitUtils::BiasType, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, MuScleFitUtils::deltaPhi(), MuScleFitUtils::deltaPhiMaxCut_, MuScleFitUtils::deltaPhiMinCut_, HLT_2023v12_cff::eta1, HLT_2023v12_cff::eta2, MuScleFitPlotter::fillTreeGen(), MuScleFitPlotter::fillTreeRec(), MuScleFitBase::genMuonPairs_, MuScleFitUtils::genMuscleFitPair, MuScleFitUtils::genPair, inputRootTreeFileName_, Data_TkAlMinBias_Run2018C_PromptReco_v3_cff::maxEvents, MuScleFitUtils::maxMuonEtaFirstRange_, MuScleFitUtils::maxMuonEtaSecondRange_, MuScleFitUtils::maxMuonPt_, MuScleFitUtils::minMuonEtaFirstRange_, MuScleFitUtils::minMuonEtaSecondRange_, MuScleFitUtils::minMuonPt_, MuScleFitBase::muonPairs_, plotter, HLT_2023v12_cff::pt1, HLT_2023v12_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);
   }
 }

◆ selGlobalMuon()

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

protected

Function for onia selections.

Definition at line 1520 of file MuScleFit.cc.

References pat::Muon::globalTrack(), pat::Muon::innerTrack(), pat::Muon::muonID(), AlCaHLTBitMon_ParallelJobs::p, and submitPVResolutionJobs::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!
 }

◆ selTrackerMuon()

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

protected

Definition at line 1535 of file MuScleFit.cc.

References pat::Muon::innerTrack(), pat::Muon::muonID(), and AlCaHLTBitMon_ParallelJobs::p.

                                                    {
   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!
 }

◆ startingNewLoop()

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 698 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;
 }

◆ takeSelectedMuonType()

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 357 of file MuScleFit.cc.

References MuScleFitBase::theMuonType_, and pwdgSkimBPark_cfi::tracks.

                                                                                 {
   // 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()));
 }

Member Data Documentation

◆ compareToSimTracks_

bool MuScleFit::compareToSimTracks_

protected

Definition at line 290 of file MuScleFit.cc.

Referenced by duringFastLoop(), and MuScleFit().

◆ fastLoop

bool MuScleFit::fastLoop

protected

Definition at line 279 of file MuScleFit.cc.

Referenced by duringLoop(), endOfLoop(), and MuScleFit().

◆ genEvtInfoToken_

edm::EDGetTokenT<GenEventInfoProduct> MuScleFit::genEvtInfoToken_

protected

Definition at line 305 of file MuScleFit.cc.

Referenced by MuScleFit(), and selectMuons().

◆ genParticlesName_

std::string MuScleFit::genParticlesName_

protected

Definition at line 293 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ iev

int MuScleFit::iev

protected

Definition at line 287 of file MuScleFit.cc.

Referenced by duringFastLoop(), endOfLoop(), and startingNewLoop().

◆ ifGenPart

bool MuScleFit::ifGenPart

protected

Definition at line 267 of file MuScleFit.cc.

◆ ifHepMC

bool MuScleFit::ifHepMC

protected

Definition at line 266 of file MuScleFit.cc.

◆ inputRootTreeFileName_

std::string MuScleFit::inputRootTreeFileName_

protected

Definition at line 296 of file MuScleFit.cc.

Referenced by duringLoop(), endOfLoop(), MuScleFit(), selectMuons(), and ~MuScleFit().

◆ loopCounter

unsigned int MuScleFit::loopCounter

protected

Definition at line 277 of file MuScleFit.cc.

Referenced by duringFastLoop(), duringLoop(), endOfFastLoop(), MuScleFit(), and startingNewLoop().

◆ maxEventsFromRootTree_

int MuScleFit::maxEventsFromRootTree_

protected

Definition at line 300 of file MuScleFit.cc.

Referenced by endOfLoop(), and MuScleFit().

◆ maxLoopNumber

unsigned int MuScleFit::maxLoopNumber

protected

Definition at line 276 of file MuScleFit.cc.

Referenced by beginOfJobInConstructor(), checkParameters(), endOfLoop(), and MuScleFit().

◆ maxResMass_hwindow

double MuScleFit::maxResMass_hwindow[6]

protected

Definition at line 272 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ minResMass_hwindow

double MuScleFit::minResMass_hwindow[6]

protected

Definition at line 271 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ muonSelector_

std::unique_ptr<MuScleFitMuonSelector> MuScleFit::muonSelector_

protected

Definition at line 317 of file MuScleFit.cc.

Referenced by MuScleFit(), and selectMuons().

◆ negateTrigger_

bool MuScleFit::negateTrigger_

protected

Definition at line 310 of file MuScleFit.cc.

Referenced by duringLoop(), and MuScleFit().

◆ numberOfEwkZ

int MuScleFit::numberOfEwkZ

protected

Definition at line 264 of file MuScleFit.cc.

◆ numberOfSimMuons

int MuScleFit::numberOfSimMuons

protected

Definition at line 262 of file MuScleFit.cc.

◆ numberOfSimTracks

int MuScleFit::numberOfSimTracks

protected

Definition at line 261 of file MuScleFit.cc.

◆ numberOfSimVertices

int MuScleFit::numberOfSimVertices

protected

Definition at line 263 of file MuScleFit.cc.

◆ outputRootTreeFileName_

std::string MuScleFit::outputRootTreeFileName_

protected

Definition at line 298 of file MuScleFit.cc.

Referenced by MuScleFit(), and ~MuScleFit().

◆ PATmuons_

bool MuScleFit::PATmuons_

protected

Definition at line 292 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ plotter

MuScleFitPlotter* MuScleFit::plotter

protected

Definition at line 281 of file MuScleFit.cc.

Referenced by beginOfJobInConstructor(), endOfFastLoop(), and selectMuons().

◆ puInfoSrc_

edm::InputTag MuScleFit::puInfoSrc_

protected

Definition at line 314 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ puInfoToken_

edm::EDGetTokenT<std::vector<PileupSummaryInfo> > MuScleFit::puInfoToken_

protected

Definition at line 304 of file MuScleFit.cc.

Referenced by MuScleFit(), and selectMuons().

◆ recMu1

reco::Particle::LorentzVector MuScleFit::recMu1

protected

Definition at line 285 of file MuScleFit.cc.

Referenced by duringFastLoop(), and selectMuons().

◆ recMu2

reco::Particle::LorentzVector MuScleFit::recMu2

protected

Definition at line 285 of file MuScleFit.cc.

Referenced by duringFastLoop(), and selectMuons().

◆ recMuScleMu1

MuScleFitMuon MuScleFit::recMuScleMu1

protected

Definition at line 286 of file MuScleFit.cc.

Referenced by selectMuons().

◆ recMuScleMu2

MuScleFitMuon MuScleFit::recMuScleMu2

protected

Definition at line 286 of file MuScleFit.cc.

Referenced by selectMuons().

◆ saveAllToTree_

bool MuScleFit::saveAllToTree_

protected

Definition at line 311 of file MuScleFit.cc.

Referenced by MuScleFit(), and ~MuScleFit().

◆ simTracksCollection_

edm::InputTag MuScleFit::simTracksCollection_

protected

Definition at line 291 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ theService

MuonServiceProxy* MuScleFit::theService

protected

Definition at line 257 of file MuScleFit.cc.

◆ totalEvents_

int MuScleFit::totalEvents_

protected

Definition at line 288 of file MuScleFit.cc.

Referenced by duringLoop(), endOfLoop(), and ~MuScleFit().

◆ triggerPath_

std::vector<std::string> MuScleFit::triggerPath_

protected

Definition at line 309 of file MuScleFit.cc.

Referenced by duringLoop(), and MuScleFit().

◆ triggerResultsLabel_

std::string MuScleFit::triggerResultsLabel_

protected

Definition at line 307 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ triggerResultsProcess_

std::string MuScleFit::triggerResultsProcess_

protected

Definition at line 308 of file MuScleFit.cc.

Referenced by duringLoop(), and MuScleFit().

◆ triggerResultsToken_

edm::EDGetTokenT<edm::TriggerResults> MuScleFit::triggerResultsToken_

protected

Definition at line 302 of file MuScleFit.cc.

Referenced by duringLoop(), and MuScleFit().

◆ vertexSrc_

edm::InputTag MuScleFit::vertexSrc_

protected

Definition at line 315 of file MuScleFit.cc.

Referenced by MuScleFit().

◆ vertexToken_

edm::EDGetTokenT<reco::VertexCollection> MuScleFit::vertexToken_

protected

Definition at line 303 of file MuScleFit.cc.

Referenced by MuScleFit(), and selectMuons().

Public Member Functions

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ MuScleFit()

◆ ~MuScleFit()

Member Function Documentation

◆ applyBias()

◆ applySmearing()

◆ beginOfJobInConstructor()

◆ checkDeltaR()

◆ checkParameters()

◆ duringFastLoop()

◆ duringLoop()

◆ endOfFastLoop()

◆ endOfJob()

◆ endOfLoop()

◆ fillComparisonHistograms()

◆ fillMuonCollection()

◆ selectMuons() [1/2]

◆ selectMuons() [2/2]

◆ selGlobalMuon()

◆ selTrackerMuon()

◆ startingNewLoop()

◆ takeSelectedMuonType()

Member Data Documentation

◆ compareToSimTracks_

◆ fastLoop

◆ genEvtInfoToken_

◆ genParticlesName_

◆ iev

◆ ifGenPart

◆ ifHepMC

◆ inputRootTreeFileName_

◆ loopCounter

◆ maxEventsFromRootTree_

◆ maxLoopNumber

◆ maxResMass_hwindow

◆ minResMass_hwindow

◆ muonSelector_

◆ negateTrigger_

◆ numberOfEwkZ

◆ numberOfSimMuons

◆ numberOfSimTracks

◆ numberOfSimVertices

◆ outputRootTreeFileName_

◆ PATmuons_

◆ plotter

◆ puInfoSrc_

◆ puInfoToken_

◆ recMu1

◆ recMu2

◆ recMuScleMu1

◆ recMuScleMu2

◆ saveAllToTree_

◆ simTracksCollection_

◆ theService

◆ totalEvents_

◆ triggerPath_

◆ triggerResultsLabel_

◆ triggerResultsProcess_

◆ triggerResultsToken_

◆ vertexSrc_

◆ vertexToken_