#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 ()
Protected Member Functions
void	applyBias (reco::Particle::LorentzVector &mu, const int charge)
	Apply the bias if needed using the function in MuScleFitUtils.
void	applySmearing (reco::Particle::LorentzVector &mu)
	Apply the smearing if needed using the function in MuScleFitUtils.
bool	checkDeltaR (reco::Particle::LorentzVector &genMu, reco::Particle::LorentzVector &recMu)
	Check if two lorentzVector are near in deltaR.
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.
void	selectMuons (const int maxEvents, const TString &treeFileName)
void	selectMuons (const edm::Event &event)
bool	selGlobalMuon (const pat::Muon *aMuon)
	Function for onia selections.
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.
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_

Detailed Description

Analyzer of the Global muon tracks

Date:: 2010/10/22 17:48:07

Revision:: 1.42

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.

References MuScleFitUtils::backgroundHandler, beginOfJobInConstructor(), MuScleFitUtils::biasFunction, MuScleFitUtils::BiasType, checkParameters(), compareToSimTracks_, MuScleFitUtils::computeMinosErrors_, gather_cfg::cout, MuScleFitUtils::crossSectionHandler, MuScleFitUtils::debug, MuScleFitBase::debug_, MuScleFitUtils::debugMassResol_, MuScleFitUtils::doBackgroundFit, MuScleFitUtils::doCrossSectionFit, MuScleFitUtils::doResolFit, MuScleFitUtils::doScaleFit, cmsRelvalreport::exit, fastLoop, MuScleFitUtils::FitStrategy, genParticlesName_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), MuScleFitUtils::goodmuon, i, inputRootTreeFileName_, j, loopCounter, MuScleFitUtils::massWindowHalfWidth, maxEventsFromRootTree_, maxLoopNumber, MuScleFitUtils::maxMuonEtaFirstRange_, MuScleFitUtils::maxMuonEtaSecondRange_, MuScleFitUtils::maxMuonPt_, maxResMass_hwindow, MuScleFitUtils::minimumShapePlots_, MuScleFitUtils::minMuonEtaFirstRange_, MuScleFitUtils::minMuonEtaSecondRange_, MuScleFitUtils::minMuonPt_, minResMass_hwindow, muonSelector_, MuScleFitUtils::MuonType, MuScleFitUtils::MuonTypeForCheckMassWindow, negateTrigger_, MuScleFitUtils::normalizeLikelihoodByEventNumber_, outputRootTreeFileName_, MuScleFitUtils::parBgr, MuScleFitUtils::parBgrFix, MuScleFitUtils::parBgrOrder, MuScleFitUtils::parBias, MuScleFitUtils::parCrossSection, MuScleFitUtils::parCrossSectionFix, MuScleFitUtils::parCrossSectionOrder, MuScleFitUtils::parResol, MuScleFitUtils::parResolFix, MuScleFitUtils::parResolOrder, MuScleFitUtils::parScale, MuScleFitUtils::parScaleFix, MuScleFitUtils::parScaleOrder, MuScleFitUtils::parSmear, PATmuons_, Pi, MuScleFitUtils::rapidityBinsForZ_, MuScleFitUtils::resfind, MuScleFitUtils::ResMass, MuScleFitUtils::ResolFitType, MuScleFitUtils::resolutionFunction, MuScleFitUtils::resolutionFunctionForVec, resolutionFunctionService(), resolutionFunctionVecService(), saveAllToTree_, MuScleFitUtils::ScaleFitType, MuScleFitUtils::scaleFunction, MuScleFitUtils::scaleFunctionForVec, scaleFunctionService(), scaleFunctionVecService(), MuScleFitUtils::sherpa_, simTracksCollection_, MuScleFitUtils::smearFunction, smearFunctionService(), MuScleFitUtils::SmearType, MuScleFitUtils::speedup, mathSSE::sqrt(), MuScleFitUtils::startWithSimplex_, AlCaHLTBitMon_QueryRunRegistry::string, MuScleFitBase::theMuonLabel_, MuScleFitBase::theMuonType_, triggerPath_, triggerResultsLabel_, triggerResultsProcess_, MuScleFitUtils::useProbsFile_, and MuScleFitUtils::x.

                                                  :
  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::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 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(), dtT0Analyzer_cfg::rootFileName, AlCaHLTBitMon_QueryRunRegistry::string, MuScleFitBase::theFiles_, MuScleFitBase::theGenInfoRootFileName_, MuScleFitBase::theRootFileName_, and MuScleFitUtils::useProbsFile_.

Referenced by MuScleFit().

{
  if (debug_>0) std::cout << "[MuScleFit]: beginOfJob" << std::endl;
  //if(maxLoopNumber>1)
  if( MuScleFitUtils::useProbsFile_ ) {
    readProbabilityDistributionsFromFile();
  }

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

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

  std::cout << "creating plotter" << std::endl;
  plotter = new MuScleFitPlotter(theGenInfoRootFileName_);
  plotter->debug = debug_;
}

bool MuScleFit::checkDeltaR	(	reco::Particle::LorentzVector &	genMu,
		reco::Particle::LorentzVector &	recMu
	)		`[protected]`

Check if two lorentzVector are near in deltaR.

Definition at line 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_, patRefSel_triggerSelection_cff::triggerResults, triggerResultsLabel_, and triggerResultsProcess_.

{
  edm::Handle<edm::TriggerResults> triggerResults;
  event.getByLabel(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()), triggerResults);
  //event.getByLabel(InputTag(triggerResultsLabel_),triggerResults);
  bool isFired = false;

  if(triggerPath_ == "")
    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_, mergeVDriftHistosByStation::name, funct::sin(), and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by duringFastLoop().

{
  std::string name(inputName + "VSMu");
  mapHisto_["hResolPt"+name]->Fill(recMu, (-genMu.Pt()+recMu.Pt())/genMu.Pt(), charge);
  mapHisto_["hResolTheta"+name]->Fill(recMu, (-genMu.Theta()+recMu.Theta()), charge);
  mapHisto_["hResolCotgTheta"+name]->Fill(recMu,(-cos(genMu.Theta())/sin(genMu.Theta())
                                                 +cos(recMu.Theta())/sin(recMu.Theta())), charge);
  mapHisto_["hResolEta"+name]->Fill(recMu, (-genMu.Eta()+recMu.Eta()),charge);
  mapHisto_["hResolPhi"+name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);

  // Fill only if it was matched to a genMu and this muon is valid
  if( (genMu.Pt() != 0) && (recMu.Pt() != 0) ) {
    mapHisto_["hPtRecoVsPt"+inputName]->Fill(genMu.Pt(), recMu.Pt());
  }
}

template<typename T >

std::vector< 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, RPCpg::mu, metsig::muon, patZpeak::muons, and mathSSE::sqrt().

{
  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_, patZpeak::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(), AlCaHLTBitMon_ParallelJobs::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(), AlCaHLTBitMon_ParallelJobs::p, and reco::HitPattern::pixelLayersWithMeasurement().

                                                   {
  
  reco::TrackRef iTrack = aMuon->innerTrack();
  const reco::HitPattern& p = iTrack->hitPattern();

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

void MuScleFit::startingNewLoop ( unsigned int int ) [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_, loopCounter, MuScleFitUtils::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()));
}