From 14304498bb1480acb72a41840c095e9d71b7b901 Mon Sep 17 00:00:00 2001 From: Ernesto Migliore Date: Thu, 30 Jul 2015 15:02:10 +0200 Subject: [PATCH] changes to deal with the post-pythia6 GEN info --- .../bin/TreeFromDump.cc | 2 +- .../bin/ZntupleToTreeConverter.cc | 2 +- .../interface/Event.h | 6 ++- .../interface/MuonPair.h | 4 +- .../plugins/BuildFile.xml | 1 + .../plugins/MuScleFit.cc | 13 +++-- .../plugins/MuScleFitMuonSelector.cc | 53 ++++++++++++++----- .../plugins/MuScleFitMuonSelector.h | 5 +- 8 files changed, 65 insertions(+), 21 deletions(-) diff --git a/MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc b/MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc index ca19614135022..208f005637867 100644 --- a/MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc +++ b/MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc @@ -73,7 +73,7 @@ int main(int argc, char* argv[]) ss >> value[i]; // std::cout << "value["<> genValue[i]; diff --git a/MuonAnalysis/MomentumScaleCalibration/bin/ZntupleToTreeConverter.cc b/MuonAnalysis/MomentumScaleCalibration/bin/ZntupleToTreeConverter.cc index b1681aa09791b..53c087a4256cf 100644 --- a/MuonAnalysis/MomentumScaleCalibration/bin/ZntupleToTreeConverter.cc +++ b/MuonAnalysis/MomentumScaleCalibration/bin/ZntupleToTreeConverter.cc @@ -142,7 +142,7 @@ int main(int argc, char* argv[]) double muon2[3] = {(*muon2pt)[i], (*muon2eta)[i], (*muon2phi)[i]}; // pairVector.push_back( std::make_pair( fromPtEtaPhiToPxPyPz(muon1), fromPtEtaPhiToPxPyPz(muon2) ) ); - pairVector.push_back( MuonPair(fromPtEtaPhiToPxPyPz(muon1), fromPtEtaPhiToPxPyPz(muon2), MuScleFitEvent(0,0,0,0,0)) ); + pairVector.push_back( MuonPair(fromPtEtaPhiToPxPyPz(muon1), fromPtEtaPhiToPxPyPz(muon2), MuScleFitEvent(0,0,0,0,0,0)) ); } } size_t namePos = fileName.find_last_of("/"); diff --git a/MuonAnalysis/MomentumScaleCalibration/interface/Event.h b/MuonAnalysis/MomentumScaleCalibration/interface/Event.h index 4f5aa8d8d29bc..2bf1111567b31 100644 --- a/MuonAnalysis/MomentumScaleCalibration/interface/Event.h +++ b/MuonAnalysis/MomentumScaleCalibration/interface/Event.h @@ -9,14 +9,16 @@ class MuScleFitEvent : public TObject MuScleFitEvent() : fRun(0), fEvent(0), + fMCWeight(0), fTrueNumPUvtx(0), fTrueNumInteractions(0), fNpv(0) {} - MuScleFitEvent(const unsigned int initRun, const unsigned long long initEvent, const int initNPUvtx, const float initTrueNI, const int initNpv) : + MuScleFitEvent(const unsigned int initRun, const unsigned long long initEvent, const double initMCWeight, const int initNPUvtx, const float initTrueNI, const int initNpv) : fRun(initRun), fEvent(initEvent), + fMCWeight(initMCWeight), fTrueNumPUvtx(initNPUvtx), fTrueNumInteractions(initTrueNI), fNpv(initNpv) @@ -25,6 +27,7 @@ class MuScleFitEvent : public TObject // Getters UInt_t run() const {return fRun;} ULong64_t event() const {return fEvent;} + Double_t MCweight() const {return fMCWeight;} Int_t nPUvtx() const {return fTrueNumPUvtx;} Float_t nTrueInteractions() const {return fTrueNumInteractions;} UInt_t npv() const {return fNpv;} @@ -32,6 +35,7 @@ class MuScleFitEvent : public TObject UInt_t fRun; ULong64_t fEvent; + Double_t fMCWeight; Int_t fTrueNumPUvtx; Float_t fTrueNumInteractions; UInt_t fNpv; diff --git a/MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h b/MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h index 9b17c0a555667..02bb6c65940f9 100644 --- a/MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h +++ b/MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h @@ -17,7 +17,7 @@ class MuonPair : public TObject //initialize 2 object of class muon mu1(lorentzVector(0,0,0,0),-1), mu2(lorentzVector(0,0,0,0),1), - event(0,0,0,0,0) + event(0,0,0,0,0,0) {} MuonPair(const MuScleFitMuon & initMu1, const MuScleFitMuon & initMu2, const MuScleFitEvent & initMuEvt) : @@ -39,7 +39,7 @@ class MuonPair : public TObject MuScleFitMuon mu2; MuScleFitEvent event; - ClassDef(MuonPair, 3) + ClassDef(MuonPair, 4) }; ClassImp(MuonPair) diff --git a/MuonAnalysis/MomentumScaleCalibration/plugins/BuildFile.xml b/MuonAnalysis/MomentumScaleCalibration/plugins/BuildFile.xml index 87db754f04c03..bb738724de8d0 100644 --- a/MuonAnalysis/MomentumScaleCalibration/plugins/BuildFile.xml +++ b/MuonAnalysis/MomentumScaleCalibration/plugins/BuildFile.xml @@ -15,6 +15,7 @@ + diff --git a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc index 5b068ba710d50..3cc9655520971 100644 --- a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc +++ b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc @@ -155,9 +155,9 @@ #include "SimDataFormats/Vertex/interface/SimVertexContainer.h" #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" - #include "DataFormats/VertexReco/interface/Vertex.h" #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" #include "TFile.h" #include "TTree.h" @@ -940,9 +940,16 @@ void MuScleFit::selectMuons(const edm::Event & event) } } + // get the MC event weight + edm::Handle genEvtInfo; + event.getByLabel("generator", genEvtInfo); + double the_genEvtweight = 1.; + if ( genEvtInfo.isValid() ) { + the_genEvtweight = genEvtInfo->weight(); + } muonPairs_.push_back(MuonPair(MuScleFitUtils::SavedPairMuScleFitMuons.back().first, MuScleFitUtils::SavedPairMuScleFitMuons.back().second, - MuScleFitEvent(event.run(), event.id().event(), the_numPUvtx, the_TrueNumInteractions, the_NVtx) + MuScleFitEvent(event.run(), event.id().event(), the_genEvtweight, the_numPUvtx, the_TrueNumInteractions, the_NVtx) )); // Fill the internal genPair tree from the external one if( MuScleFitUtils::speedup == false ) { @@ -1036,7 +1043,7 @@ void MuScleFit::selectMuons(const int maxEvents, const TString & treeFileName) //FIXME: we loose the additional information besides the 4-momenta muonPairs_.push_back(MuonPair(MuScleFitMuon(it->first,-1), MuScleFitMuon(it->second,+1), - MuScleFitEvent((*evtRunIt).first, (*evtRunIt).second, 0, 0, 0) ) // FIXME: order of event and run number mixed up! + MuScleFitEvent((*evtRunIt).first, (*evtRunIt).second, 0, 0, 0, 0) ) // FIXME: order of event and run number mixed up! ); diff --git a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.cc b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.cc index 3439d2b3cbd70..d73965fb1e97c 100644 --- a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.cc +++ b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.cc @@ -8,23 +8,44 @@ const unsigned int MuScleFitMuonSelector::motherPdgIdArray[] = {23, 100553, 1005 const reco::Candidate* MuScleFitMuonSelector::getStatus1Muon(const reco::Candidate* status3Muon){ - const reco::Candidate* tempStatus1Muon = status3Muon; - int status = tempStatus1Muon->status(); - while(tempStatus1Muon == 0 || tempStatus1Muon->numberOfDaughters()!=0){ - if (status == 1) break; + const reco::Candidate* tempMuon = status3Muon; + // bool lastCopy = ((reco::GenParticle*)tempMuon)->isLastCopy(); // isLastCopy() likely not enough robust + bool isPromptFinalState = ((reco::GenParticle*)tempMuon)->isPromptFinalState(); // pre-CMSSW_74X code: int status = tempStatus1Muon->status(); + while(tempMuon == 0 || tempMuon->numberOfDaughters()!=0){ + if ( isPromptFinalState ) break; // pre-CMSSW_74X code: if (status == 1) break; //std::vector daughters; - for (unsigned int i=0; inumberOfDaughters(); ++i){ - if ( tempStatus1Muon->daughter(i)->pdgId()==tempStatus1Muon->pdgId() ){ - tempStatus1Muon = tempStatus1Muon->daughter(i); - status = tempStatus1Muon->status(); + for (unsigned int i=0; inumberOfDaughters(); ++i){ + if ( tempMuon->daughter(i)->pdgId()==tempMuon->pdgId() ){ + tempMuon = tempMuon->daughter(i); + isPromptFinalState = ((reco::GenParticle*)tempMuon)->isPromptFinalState(); // pre-CMSSW_74X code: status = tempStatus1Muon->status(); break; }else continue; }//for loop }//while loop - return tempStatus1Muon; + return tempMuon; } +const reco::Candidate* +MuScleFitMuonSelector::getStatus3Muon(const reco::Candidate* status3Muon){ + const reco::Candidate* tempMuon = status3Muon; + bool lastCopy = ((reco::GenParticle*)tempMuon)->isLastCopyBeforeFSR(); // pre-CMSSW_74X code: int status = tempStatus1Muon->status(); + while(tempMuon == 0 || tempMuon->numberOfDaughters()!=0){ + if ( lastCopy ) break; // pre-CMSSW_74X code: if (status == 3) break; + //std::vector daughters; + for (unsigned int i=0; inumberOfDaughters(); ++i){ + if ( tempMuon->daughter(i)->pdgId()==tempMuon->pdgId() ){ + tempMuon = tempMuon->daughter(i); + lastCopy = ((reco::GenParticle*)tempMuon)->isLastCopyBeforeFSR(); // pre-CMSSW_74X code: status = tempStatus1Muon->status(); + break; + }else continue; + }//for loop + }//while loop + + return tempMuon; +} + + bool MuScleFitMuonSelector::selGlobalMuon(const pat::Muon* aMuon) { @@ -386,7 +407,7 @@ GenMuonPair MuScleFitMuonSelector::findGenMuFromRes( const edm::HepMCProduct* ev muFromRes.mu2 = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(), (*part)->momentum().pz(),(*part)->momentum().e())); } - muFromRes.motherId = motherPdgId; + muFromRes.motherId = motherPdgId; } } } @@ -410,7 +431,14 @@ GenMuonPair MuScleFitMuonSelector::findGenMuFromRes( const reco::GenParticleColl std::cout << "Found a muon with mother: " << motherPdgId << std::endl; } for( int ires = 0; ires < 6; ++ires ) { - if( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) fromRes = true; + // if( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) fromRes = true; // changed by EM 2015.07.30 + // begin of comment + // the list of resonances motherPdgIdArray does not contain the photon (PdgId = 21) while ~1% of the + // mu+mu- events in the range [50,120] GeV has a photon as the mother. + // It needs to be fixed without spoiling the logic of the selection of different resonances + // e.g. mixing up onia etc. + // end of comment + if( ( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) || motherPdgId == 21 ) fromRes = true; } if(fromRes){ if (debug_>0) std::cout<<"fromRes = true, motherPdgId = "<p4().px(),status1Muon->p4().py(), // status1Muon->p4().pz(),status1Muon->p4().e())); } + muFromRes.motherId = motherPdgId; } } } @@ -455,7 +484,7 @@ std::pair MuScleFitMuonSelector::findSimMuFromRes( bool fromRes = false; unsigned int motherPdgId = (*mother)->pdg_id(); for( int ires = 0; ires < 6; ++ires ) { - if( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) fromRes = true; + if( ( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) || motherPdgId == 21 ) fromRes = true; } if( fromRes ) { if(gp->pdg_id() == 13) diff --git a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.h b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.h index 8ed68bbca3ad5..ff80a2fab531c 100644 --- a/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.h +++ b/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.h @@ -47,8 +47,11 @@ class MuScleFitMuonSelector {} ~MuScleFitMuonSelector() {} - //Method to get the muon after FSR (status 1 muon) starting from status 3 muon which is daughter of the Z + //Method to get the muon after FSR (status 1 muon in PYTHIA6) starting from status 3 muon which is daughter of the Z const reco::Candidate* getStatus1Muon(const reco::Candidate* status3Muon); + + //Method to get the muon before FSR (status 3 muon in PYTHIA6) starting from status 3 muon which is daughter of the Z + const reco::Candidate* getStatus3Muon(const reco::Candidate* status3Muon); /// Main method used to select muons of type specified by muonType_ from the collection specified by muonLabel_ and PATmuons_ void selectMuons(const edm::Event & event, std::vector & muons,