#include <MuScleFit.h>

Inheritance diagram for MuScleFit:

Public Member Functions
void	beginOfJobInConstructor ()

virtual void	duringFastLoop ()

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

virtual void	endOfFastLoop (const unsigned int iLoop)

virtual void	endOfJob ()

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

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

	MuScleFit (const edm::ParameterSet &pset)

virtual void	startingNewLoop (unsigned int iLoop)

virtual	~MuScleFit ()

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

virtual	~EDLooper ()

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...

virtual void	beginOfJob (EventSetup const &)

virtual void	beginOfJob ()

void	copyInfo (ScheduleInfo const &)

void	doBeginLuminosityBlock (LuminosityBlockPrincipal &, EventSetup const &)

void	doBeginRun (RunPrincipal &, EventSetup const &)

Status	doDuringLoop (EventPrincipal &eventPrincipal, EventSetup const &es, ProcessingController &)

void	doEndLuminosityBlock (LuminosityBlockPrincipal &, EventSetup const &)

Status	doEndOfLoop (EventSetup const &es)

void	doEndRun (RunPrincipal &, EventSetup const &)

void	doStartingNewLoop ()

	EDLooperBase ()

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

void	prepareForNextLoop (eventsetup::EventSetupProvider *esp)

void	setActionTable (ActionTable const *actionTable)

void	setModuleChanger (ModuleChanger const *)

virtual	~EDLooperBase ()

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 const *	moduleChanger () const
	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::auto_ptr < MuScleFitMuonSelector >	muonSelector_

bool	negateTrigger_

int	numberOfEwkZ

int	numberOfSimMuons

int	numberOfSimTracks

int	numberOfSimVertices

std::string	outputRootTreeFileName_

bool	PATmuons_

MuScleFitPlotter *	plotter

reco::Particle::LorentzVector	recMu1

reco::Particle::LorentzVector	recMu2

bool	saveAllToTree_

edm::InputTag	simTracksCollection_

MuonServiceProxy *	theService

int	totalEvents_

std::string	triggerPath_

std::string	triggerResultsLabel_

std::string	triggerResultsProcess_

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 ()

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

Date:: 2010/09/28 14:45:38

Revision:: 1.40

Author: C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo - INFN Padova

Definition at line 50 of file MuScleFit.h.

Constructor & Destructor Documentation

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

Definition at line 148 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] = 76.;
   maxResMass_hwindow[0] = 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::parScale        = pset.getParameter<std::vector<double> >("parScale");
   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::string>("TriggerPath");
   negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
   saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
 
   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::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::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_ == 4 || 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.5;
   MuScleFitUtils::massWindowHalfWidth[0][2] = 0.5;
   MuScleFitUtils::massWindowHalfWidth[0][3] = 0.5;
   MuScleFitUtils::massWindowHalfWidth[0][4] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[0][5] = 0.2;
   MuScleFitUtils::massWindowHalfWidth[1][0] = 50.;
   MuScleFitUtils::massWindowHalfWidth[1][1] = 2.5;
   MuScleFitUtils::massWindowHalfWidth[1][2] = 2.5;
   MuScleFitUtils::massWindowHalfWidth[1][3] = 2.5;
   MuScleFitUtils::massWindowHalfWidth[1][4] = 1.5;
   MuScleFitUtils::massWindowHalfWidth[1][5] = 1.5;
   MuScleFitUtils::massWindowHalfWidth[2][0] = 20.;
   MuScleFitUtils::massWindowHalfWidth[2][1] = 0.5;
   MuScleFitUtils::massWindowHalfWidth[2][2] = 0.5;
   MuScleFitUtils::massWindowHalfWidth[2][3] = 0.5;
   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 ( )

virtual

Definition at line 357 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, MuScleFitBase::genMuonPairs_, inputRootTreeFileName_, 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_);
         rootTreeHandler.writeTree(outputRootTreeFileName_, &(muonPairs_), theMuonType_, 0, 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 948 of file MuScleFit.cc.

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

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

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

protected

Apply the smearing if needed using the function in MuScleFitUtils.

Definition at line 939 of file MuScleFit.cc.

References MuScleFitUtils::applySmearing(), gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::goodmuon, 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;
   }
 }

void MuScleFit::beginOfJobInConstructor ( )

Definition at line 383 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitPlotter::debug, MuScleFitBase::debug_, i, maxLoopNumber, plotter, MuScleFitBase::readProbabilityDistributionsFromFile(), 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 908 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::deltaPhi(), 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 959 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.size()!=0) ||
       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 682 of file MuScleFit.cc.

References CrossSectionHandler::addParameters(), MuScleFitUtils::applyScale(), checkDeltaR(), compareToSimTracks_, MuScleFitUtils::computeWeight(), gather_cfg::cout, MuScleFitUtils::crossSectionHandler, MuScleFitBase::debug_, MuScleFitUtils::debugMassResol_, MuScleFitUtils::doScaleFit, fillComparisonHistograms(), first, MuScleFitUtils::genPair, i, iev, MuScleFitUtils::iev_, funct::log(), loopCounter, MuScleFitBase::mapHisto_, MuScleFitUtils::massProb(), MuScleFitUtils::massResolComponents, MuScleFitUtils::massResolution(), MuScleFitUtils::parBgr, MuScleFitUtils::parCrossSection, MuScleFitUtils::parResol, MuScleFitUtils::parScale, MuScleFitUtils::parvalue, funct::pow(), recMu1, recMu2, MuScleFitUtils::ResFound, MuScleFitUtils::resolutionFunctionForVec, MuScleFitUtils::SavedPair, edm::second(), MuScleFitUtils::simPair, findQualityFiles::size, MuScleFitUtils::speedup, funct::true, and CommonMethods::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);
   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);
     mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
     mapHisto_["hRecBestMu"]->Fill(recMu2);
     mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
     mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
     // Reconstructed resonance
     mapHisto_["hRecBestRes"]->Fill(bestRecRes, weight);
     mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, 1.);
     // Fill histogram of Res mass vs muon variables
     mapHisto_["hRecBestResVSMu"]->Fill (recMu1, bestRecRes, -1);
     mapHisto_["hRecBestResVSMu"]->Fill (recMu2, bestRecRes, +1);
     // Fill histogram of Res mass vs Res variables
     mapHisto_["hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
 
     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 );
       } 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 );
       }
       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 &
	)

virtual

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

Implements edm::EDLooper.

Definition at line 528 of file MuScleFit.cc.

References gather_cfg::cout, MuScleFitBase::debug_, duringFastLoop(), fastLoop, edm::Event::getRun(), HltComparatorCreateWorkflow::hltConfig, HLTConfigProvider::init(), inputRootTreeFileName_, edm::EDLooperBase::kContinue, loopCounter, negateTrigger_, selectMuons(), findQualityFiles::size, totalEvents_, HLTConfigProvider::triggerIndex(), triggerPath_, 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_ == "")
     isFired = true;
   else if(triggerPath_ == "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);
     unsigned int triggerIndex( hltConfig.triggerIndex(triggerPath_) );
     // 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_<<" "<<isFired<<std::endl;
     }
   }
 
   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 492 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 ( )

virtual

Reimplemented from edm::EDLooperBase.

Definition at line 412 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 &	,
		unsigned int	iCounter
	)

virtual

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 445 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%1000 == 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 922 of file MuScleFit.cc.

References DeDxDiscriminatorTools::charge(), funct::cos(), MuScleFitUtils::deltaPhiNoFabs(), MuScleFitBase::mapHisto_, AlCaRecoCosmics_cfg::name, and funct::sin().

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< reco::LeafCandidate > MuScleFit::fillMuonCollection ( const std::vector< T > & tracks )

Definition at line 174 of file MuScleFit.h.

References applyBias(), applySmearing(), gather_cfg::cout, MuScleFitBase::debug_, MuScleFitUtils::goodmuon, MuScleFitUtils::mMu2, metsig::muon, ExpressReco_HICollisions_FallBack::muons, mathSSE::sqrt(), and ExpressReco_HICollisions_FallBack::track.

 {
   std::vector<reco::LeafCandidate> 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());
 
     reco::LeafCandidate muon(track->charge(),mu);
     // 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 590 of file MuScleFit.cc.

References prof2calltree::back, gather_cfg::cout, MuScleFitBase::debug_, MuScleFitPlotter::fillRec(), MuScleFitUtils::findBestRecoRes(), first, MuScleFitBase::genMuonPairs_, MuScleFitUtils::genPair, MuScleFitBase::muonPairs_, ExpressReco_HICollisions_FallBack::muons, muonSelector_, plotter, reco::tau::disc::Pt(), recMu1, recMu2, MuScleFitUtils::ResFound, edm::Event::run(), MuScleFitUtils::SavedPair, edm::second(), and MuScleFitUtils::simPair.

Referenced by duringLoop(), and endOfLoop().

 {
   recMu1 = reco::Particle::LorentzVector(0,0,0,0);
   recMu2 = reco::Particle::LorentzVector(0,0,0,0);
 
   std::vector<reco::LeafCandidate> muons;
   muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
   plotter->fillRec(muons);
 
   // Find the two muons from the resonance, and set ResFound bool
   // ------------------------------------------------------------
   std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> 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;
     recMu2 = 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;
     }
     MuScleFitUtils::SavedPair.push_back( std::make_pair( recMu1, recMu2 ) );
   } else {
     MuScleFitUtils::SavedPair.push_back( std::make_pair( lorentzVector(0.,0.,0.,0.), lorentzVector(0.,0.,0.,0.) ) );
   }
   // Save the events also in the external tree so that it can be saved later
   muonPairs_.push_back(MuonPair(MuScleFitUtils::SavedPair.back().first,
                 MuScleFitUtils::SavedPair.back().second,
                 event.run(), event.id().event()));
   // Fill the internal genPair tree from the external one
   MuScleFitUtils::genPair.push_back(std::make_pair( genMuonPairs_.back().mu1, genMuonPairs_.back().mu2 ));
 }

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

References applyBias(), applySmearing(), MuScleFitUtils::BiasType, gather_cfg::cout, MuScleFitPlotter::fillGen(), MuScleFitPlotter::fillRec(), MuScleFitUtils::genPair, inputRootTreeFileName_, MuScleFitUtils::maxMuonEtaFirstRange_, MuScleFitUtils::maxMuonEtaSecondRange_, MuScleFitUtils::maxMuonPt_, MuScleFitUtils::minMuonEtaFirstRange_, MuScleFitUtils::minMuonEtaSecondRange_, MuScleFitUtils::minMuonPt_, plotter, RootTreeHandler::readTree(), MuScleFitUtils::SavedPair, MuScleFitUtils::SmearType, MuScleFitUtils::speedup, and MuScleFitBase::theMuonType_.

 {
   std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
   RootTreeHandler rootTreeHandler;
   if( MuScleFitUtils::speedup ) {
     rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_);
   }
   else {
     rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, &(MuScleFitUtils::genPair));
   }
   // Now loop on all the pairs and apply any smearing and bias if needed
   std::vector<std::pair<lorentzVector,lorentzVector> >::iterator it = MuScleFitUtils::SavedPair.begin();
   for( ; it != MuScleFitUtils::SavedPair.end(); ++it ) {
 
     // 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
     if( !(pt1 > MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
       pt2 > MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
       ( (eta1 > MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_ &&
          eta2 > MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_) ||
         (eta1 > MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_ &&
          eta2 > MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_) ) ) ) {
       // std::cout << "removing muons not passing cuts" << std::endl;
       it->first = reco::Particle::LorentzVector(0,0,0,0);
       it->second = reco::Particle::LorentzVector(0,0,0,0);
     }
 
     // First is always mu-, second mu+
     if( (MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0) ) {
       applySmearing(it->first);
       applyBias(it->first, -1);
       applySmearing(it->second);
       applyBias(it->second, 1);
     }
   }
   plotter->fillRec(MuScleFitUtils::SavedPair);
   if( !(MuScleFitUtils::speedup) ) {
     plotter->fillGen(MuScleFitUtils::genPair);
   }
 }

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

protected

Function for onia selections.

Definition at line 1046 of file MuScleFit.cc.

References pat::Muon::globalTrack(), pat::Muon::innerTrack(), pat::Muon::muonID(), reco::HitPattern::numberOfValidMuonHits(), L1TEmulatorMonitor_cff::p, reco::HitPattern::pixelLayersWithMeasurement(), and lumiQueryAPI::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 1067 of file MuScleFit.cc.

References pat::Muon::innerTrack(), pat::Muon::muonID(), L1TEmulatorMonitor_cff::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 )

virtual

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 418 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 201 of file MuScleFit.h.

References MuScleFitBase::theMuonType_.

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