111 #include "HepPDT/defs.h"
112 #include "HepPDT/TableBuilder.hh"
113 #include "HepPDT/ParticleDataTable.hh"
115 #include "HepMC/GenParticle.h"
116 #include "HepMC/GenEvent.h"
136 #include "valgrind/callgrind.h"
155 if ((theMuonType_<-4 || theMuonType_>5) &&
theMuonType_<10) {
156 std::cout <<
"[MuScleFit]: Unknown muon type! Aborting." << std::endl;
200 int resolFitType = pset.
getParameter<
int>(
"ResolFitType");
272 if (MuScleFitUtils::SmearType>0) {
273 std::cout <<
"[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
274 TF1 G(
"G",
"[0]*exp(-0.5*pow(x,2))", -5., 5.);
276 G.SetParameter (0,norm);
277 for (
int i=0;
i<10000;
i++) {
278 for (
int j=0;
j<7;
j++) {
300 MuScleFitUtils::rapidityBinsForZ_ =
false;
325 MuScleFitUtils::resfind,
326 MuScleFitUtils::speedup, genParticlesName_,
327 compareToSimTracks_, simTracksCollection_,
328 MuScleFitUtils::sherpa_,
debug_));
331 pset.
getParameter<std::vector<double> >(
"LeftWindowBorder"),
332 pset.
getParameter<std::vector<double> >(
"RightWindowBorder"),
364 std::cout <<
"Saving muon pairs to root tree" << std::endl;
398 std::stringstream ss;
401 theFiles_.push_back (
new TFile(rootFileName.c_str(),
"RECREATE"));
403 if (
debug_>0)
std::cout <<
"[MuScleFit]: Root file created" << std::endl;
405 std::cout <<
"creating plotter" << std::endl;
420 if (
debug_>0)
std::cout <<
"[MuScleFit]: Starting loop # " << iLoop << std::endl;
447 unsigned int iFastLoop = 1;
464 std::cout <<
"Starting fast loop number " << iFastLoop << std::endl;
474 if(
iev%1000 == 0 ) {
475 std::cout <<
"Fast looping on event number " <<
iev << std::endl;
480 std::cout <<
"End of fast loop number " << iFastLoop <<
". Ran on " <<
iev <<
" events" << std::endl;
502 std::cout <<
"Ending loop # " << iLoop << std::endl;
512 TDirectory * likelihoodDir =
theFiles_[iLoop]->mkdir(
"likelihood");
533 bool isFired =
false;
538 isFired =triggerResults->accept();
540 std::cout<<
"Trigger "<<isFired<<std::endl;
548 if (triggerIndex < triggerResults->
size()) {
549 isFired = triggerResults->accept(triggerIndex);
559 CALLGRIND_START_INSTRUMENTATION;
564 <<
" Run: " <<
event.id().run() <<
" Event: " <<
event.id().event() << std::endl;
585 CALLGRIND_STOP_INSTRUMENTATION;
586 CALLGRIND_DUMP_STATS;
595 std::vector<reco::LeafCandidate>
muons;
601 std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> recMuFromBestRes =
606 std::cout <<std::setprecision(9)<<
"Pt after findbestrecores: " << (recMuFromBestRes.first).
Pt() <<
" "
607 << (recMuFromBestRes.second).
Pt() << std::endl;
611 recMu1 = recMuFromBestRes.first;
612 recMu2 = recMuFromBestRes.second;
626 event.
run(),
event.id().event()));
633 std::cout <<
"Reading muon pairs from Root Tree in " << treeFileName << std::endl;
648 double pt1 = it->first.pt();
650 double pt2 = it->second.pt();
652 double eta1 = it->first.eta();
654 double eta2 = it->second.eta();
695 if(
debug_>0 )
std::cout <<
"About to start lik par correction and histo filling; ResFound is "
708 <<
recMu1.Pt() <<
" Pt2 = " <<
recMu2.Pt() << std::endl;
720 <<
recMu1.Pt() <<
" Pt2 = " <<
recMu2.Pt() << std::endl;
733 mapHisto_[
"hRecBestRes"]->Fill(bestRecRes, weight);
734 mapHisto_[
"hRecBestResAllEvents"]->Fill(bestRecRes, 1.);
739 mapHisto_[
"hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
741 std::vector<double> * parval;
742 std::vector<double> initpar;
790 if(
debug_ > 0 )
std::cout <<
"mass = " << bestRecRes.mass() << std::endl;
796 std::vector<double> initpar;
819 if(
debug_ > 0 )
std::cout <<
"inside weight: mass = " << bestRecRes.mass() <<
", prob = " << prob << std::endl;
821 if(
debug_ > 0 )
std::cout <<
"inside prob: mass = " << bestRecRes.mass() <<
", prob = " << prob << std::endl;
830 if( recoMass != 0 ) {
834 mapHisto_[
"hFunctionResolMassVSMu"]->Fill(
recMu1, massResol/recoMass, -1);
835 mapHisto_[
"hFunctionResolMassVSMu"]->Fill(
recMu2, massResol/recoMass, +1);
853 double diffMass = (recoMass - genMass)/genMass;
857 double eta1 =
recMu1.eta();
859 double eta2 =
recMu2.eta();
861 if( diffMass == diffMass ) {
863 mapHisto_[
"hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
864 mapHisto_[
"hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
865 mapHisto_[
"hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
866 mapHisto_[
"hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
868 mapHisto_[
"hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
869 mapHisto_[
"hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
872 std::cout <<
"Error, there is a nan: recoMass = " << recoMass <<
", genMass = " << genMass << std::endl;
881 mapHisto_[
"hMass_P"]->Fill(bestRecRes.mass(), prob);
882 if(
debug_ > 0 )
std::cout <<
"mass = " << bestRecRes.mass() <<
", prob = " << prob << std::endl;
883 mapHisto_[
"hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
885 mapHisto_[
"hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
888 mapHisto_[
"hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
898 if (
debug_>0)
std::cout <<
"[MuScleFit]: filling the pair" << std::endl;
911 ((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
915 std::cout<<
"Reco muon "<<recMu<<
" with eta "<<recMu.Eta()<<
" and phi "<<recMu.Phi()<<std::endl
916 <<
" DOES NOT MATCH with generated muon from resonance: "<<std::endl
917 <<genMu<<
" with eta "<<genMu.Eta()<<
" and phi "<<genMu.Phi()<<std::endl;
923 const std::string & inputName,
const int charge )
925 std::string
name(inputName +
"VSMu");
927 mapHisto_[
"hResolTheta"+
name]->Fill(recMu, (-genMu.Theta()+recMu.Theta()), charge);
930 mapHisto_[
"hResolEta"+
name]->Fill(recMu, (-genMu.Eta()+recMu.Eta()),charge);
934 if( (genMu.Pt() != 0) && (recMu.Pt() != 0) ) {
935 mapHisto_[
"hPtRecoVsPt"+inputName]->Fill(genMu.Pt(), recMu.Pt());
944 <<
": after smearing Pt = " << mu.Pt() << std::endl;
953 <<
": after bias Pt = " << mu.Pt() << std::endl;
964 std::cout <<
"it must have as many values as the number of loops, which is = " <<
maxLoopNumber << std::endl;
969 std::cout <<
"it must have as many values as the number of loops, which is = " <<
maxLoopNumber << std::endl;
974 std::cout <<
"it must have as many values as the number of loops, which is = " <<
maxLoopNumber << std::endl;
979 std::cout <<
"it must have as many values as the number of loops, which is = " <<
maxLoopNumber << std::endl;
998 MuScleFitUtils::BiasType<0 || MuScleFitUtils::BiasType>13) {
999 std::cout <<
"[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
1011 MuScleFitUtils::SmearType<0 || MuScleFitUtils::SmearType>7) {
1012 std::cout <<
"[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
1019 std::cout <<
"[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
1024 std::cout <<
"[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
1029 std::cout <<
"[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1034 std::cout <<
"[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1041 std::cout <<
"[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
1055 iTrack->found() > 11 &&
1056 gTrack->chi2()/gTrack->ndof() < 20.0 &&
1058 iTrack->chi2()/iTrack->ndof() < 4.0 &&
1059 aMuon->
muonID(
"TrackerMuonArbitrated") &&
1060 aMuon->
muonID(
"TMLastStationAngTight") &&
1062 fabs(iTrack->dxy()) < 3.0 &&
1063 fabs(iTrack->dz()) < 15.0 );
1073 iTrack->found() > 11 &&
1074 iTrack->chi2()/iTrack->ndof() < 4.0 &&
1075 aMuon->
muonID(
"TrackerMuonArbitrated") &&
1076 aMuon->
muonID(
"TMLastStationAngTight") &&
1078 fabs(iTrack->dxy()) < 3.0 &&
1079 fabs(iTrack->dz()) < 15.0 );
static double deltaPhiNoFabs(const double &phi1, const double &phi2)
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance c...
static std::vector< std::pair< lorentzVector, lorentzVector > > simPair
static std::vector< int > doScaleFit
static std::vector< int > doResolFit
T getParameter(std::string const &) const
scaleFunctionBase< double * > * scaleFunctionService(const int identifier)
Service to build the scale functor corresponding to the passed identifier.
T getUntrackedParameter(std::string const &, T const &) const
static std::vector< double > parBias
void selectMuons(const edm::Event &event)
virtual edm::EDLooper::Status endOfLoop(const edm::EventSetup &eventSetup, unsigned int iLoop)
static std::vector< int > parCrossSectionOrder
std::vector< GenMuonPair > genMuonPairs_
Stores the genMuon pairs and the motherId prior to the creation of the internal tree.
reco::TrackRef innerTrack() const
reference to Track reconstructed in the tracker only (reimplemented from reco::Muon) ...
static smearFunctionBase * smearFunction
static std::vector< int > parScaleOrder
edm::InputTag theMuonLabel_
static std::vector< int > doCrossSectionFit
static void minimizeLikelihood()
static double maxMuonEtaSecondRange_
static std::vector< int > parBgrOrder
bool selGlobalMuon(const pat::Muon *aMuon)
Function for onia selections.
static unsigned int loopCounter
bool muonID(const std::string &name) const
int numberOfValidMuonHits() const
static bool startWithSimplex_
std::string triggerResultsProcess_
static std::vector< int > doBackgroundFit
void fillRec(std::vector< reco::LeafCandidate > &muons)
static std::vector< double > parResol
static bool debugMassResol_
virtual void duringFastLoop()
Sin< T >::type sin(const T &t)
static BackgroundHandler * backgroundHandler
static double x[7][10000]
static double massWindowHalfWidth[3][6]
Run const & getRun() const
std::map< std::string, Histograms * > mapHisto_
The map of histograms.
static scaleFunctionBase< std::vector< double > > * biasFunction
bool selTrackerMuon(const pat::Muon *aMuon)
static std::vector< int > parBgrFix
static bool minimumShapePlots_
int pixelLayersWithMeasurement() const
edm::InputTag simTracksCollection_
static int MuonTypeForCheckMassWindow
static double minMuonEtaFirstRange_
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
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.
static struct MuScleFitUtils::massResolComponentsStruct massResolComponents
static std::vector< int > parScaleFix
static bool computeMinosErrors_
reco::Particle::LorentzVector lorentzVector
static scaleFunctionBase< std::vector< double > > * scaleFunctionForVec
static std::vector< std::vector< double > > parvalue
unsigned int triggerIndex(const std::string &triggerName) const
slot position of trigger path in trigger table (0 to size-1)
U second(std::pair< T, U > const &p)
static std::vector< std::pair< lorentzVector, lorentzVector > > genPair
std::string outputRootTreeFileName_
std::string theGenInfoRootFileName_
resolutionFunctionBase< std::vector< double > > * resolutionFunctionVecService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier when receiving a std::...
virtual void startingNewLoop(unsigned int iLoop)
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
reco::Particle::LorentzVector recMu2
reco::Particle::LorentzVector recMu1
static lorentzVector applyBias(const lorentzVector &muon, const int charge)
static std::vector< double > parBgr
unsigned int maxLoopNumber
std::string genParticlesName_
Cos< T >::type cos(const T &t)
void clearHistoMap()
Clean the histograms map.
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, MuonPairVector *genPair=0)
static double computeWeight(const double &mass, const int iev, const bool doUseBkgrWindow=false)
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
MuScleFit(const edm::ParameterSet &pset)
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)
static std::vector< int > parResolFix
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
static std::vector< double > parSmear
void beginOfJobInConstructor()
smearFunctionBase * smearFunctionService(const int identifier)
Service to build the smearing functor corresponding to the passed identifier.
reco::TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector (reimplemented from reco::Muon) ...
static std::vector< int > parResolOrder
std::auto_ptr< MuScleFitMuonSelector > muonSelector_
double deltaR(double eta1, double eta2, double phi1, double phi2)
void writeTree(const TString &fileName, const std::vector< MuonPair > *savedPair, const int muonType=0, const std::vector< GenMuonPair > *genPair=0, const bool saveAll=false)
scaleFunctionBase< std::vector< double > > * scaleFunctionVecService(const int identifier)
Service to build the scale functor corresponding to the passed identifier when receiving a std::vecto...
void writeHistoMap(const unsigned int iLoop)
Save the histograms map to file.
static std::vector< double > parScale
static double oldNormalization_
void applySmearing(reco::Particle::LorentzVector &mu)
Apply the smearing if needed using the function in MuScleFitUtils.
Log< T >::type log(const T &t)
std::string theRootFileName_
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
std::string triggerResultsLabel_
static bool rapidityBinsForZ_
MuScleFitPlotter * plotter
void fillHistoMap(TFile *outputFile, unsigned int iLoop)
Create the histograms map.
std::vector< TFile * > theFiles_
The files were the histograms are saved.
double minResMass_hwindow[6]
static double minMuonEtaSecondRange_
void fillGen(const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &genPairs)
static std::vector< double > parCrossSection
static lorentzVector applySmearing(const lorentzVector &muon)
static bool normalizeLikelihoodByEventNumber_
std::string inputRootTreeFileName_
bool checkDeltaR(reco::Particle::LorentzVector &genMu, reco::Particle::LorentzVector &recMu)
Check if two lorentzVector are near in deltaR.
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
int maxEventsFromRootTree_
double maxResMass_hwindow[6]
std::vector< MuonPair > muonPairs_
Used to store the muon pairs plus run and event number prior to the creation of the internal tree...
static double deltaPhi(const double &phi1, const double &phi2)
static scaleFunctionBase< double * > * scaleFunction
static bool useProbsFile_
static std::vector< int > resfind
void readProbabilityDistributionsFromFile()
Read probability distributions from a local root file.
static int counter_resprob
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Analysis-level muon class.
static double maxMuonEtaFirstRange_
tuple size
Write out results.
static std::pair< lorentzVector, lorentzVector > findBestRecoRes(const std::vector< reco::LeafCandidate > &muons)
Power< A, B >::type pow(const A &a, const B &b)
void applyBias(reco::Particle::LorentzVector &mu, const int charge)
Apply the bias if needed using the function in MuScleFitUtils.
static CrossSectionHandler * crossSectionHandler
static std::vector< int > parCrossSectionFix
static resolutionFunctionBase< double * > * resolutionFunction
virtual edm::EDLooper::Status duringLoop(const edm::Event &event, const edm::EventSetup &eventSetup)
void addParameters(std::vector< double > &initpar)
Inputs the vars in a vector.
virtual void endOfFastLoop(const unsigned int iLoop)