Skip to content

Commit

Permalink
changes to deal with the post-pythia6 GEN info
Browse files Browse the repository at this point in the history
  • Loading branch information
emiglior committed Jul 30, 2015
1 parent 78b3b2c commit 1430449
Show file tree
Hide file tree
Showing 8 changed files with 65 additions and 21 deletions.
2 changes: 1 addition & 1 deletion MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc
Expand Up @@ -73,7 +73,7 @@ int main(int argc, char* argv[])
ss >> value[i];
// std::cout << "value["<<i<<"] = " << value[i] << std::endl;
}
pairVector.push_back(MuonPair(fromPtEtaPhiToPxPyPz(value), fromPtEtaPhiToPxPyPz(&(value[3])), MuScleFitEvent(0,0,0,0,0)) );
pairVector.push_back(MuonPair(fromPtEtaPhiToPxPyPz(value), fromPtEtaPhiToPxPyPz(&(value[3])), MuScleFitEvent(0,0,0,0,0,0)) );
if( genInfo ) {
for( int i=0; i<6; ++i ) {
ss >> genValue[i];
Expand Down
Expand Up @@ -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("/");
Expand Down
6 changes: 5 additions & 1 deletion MuonAnalysis/MomentumScaleCalibration/interface/Event.h
Expand Up @@ -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)
Expand All @@ -25,13 +27,15 @@ 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;}


UInt_t fRun;
ULong64_t fEvent;
Double_t fMCWeight;
Int_t fTrueNumPUvtx;
Float_t fTrueNumInteractions;
UInt_t fNpv;
Expand Down
4 changes: 2 additions & 2 deletions MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h
Expand Up @@ -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) :
Expand All @@ -39,7 +39,7 @@ class MuonPair : public TObject
MuScleFitMuon mu2;
MuScleFitEvent event;

ClassDef(MuonPair, 3)
ClassDef(MuonPair, 4)
};
ClassImp(MuonPair)

Expand Down
Expand Up @@ -15,6 +15,7 @@
<use name="SimGeneral/HepPDTRecord"/>
<use name="SimDataFormats/Track"/>
<use name="SimDataFormats/Vertex"/>
<use name="SimDataFormats/GeneratorProducts"/>
<use name="root"/>
<use name="clhep"/>
<use name="heppdt"/>
Expand Down
13 changes: 10 additions & 3 deletions MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc
Expand Up @@ -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"
Expand Down Expand Up @@ -940,9 +940,16 @@ void MuScleFit::selectMuons(const edm::Event & event)
}
}

// get the MC event weight
edm::Handle<GenEventInfoProduct> 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 ) {
Expand Down Expand Up @@ -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!
);


Expand Down
Expand Up @@ -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<const reco::Candidate*> daughters;
for (unsigned int i=0; i<tempStatus1Muon->numberOfDaughters(); ++i){
if ( tempStatus1Muon->daughter(i)->pdgId()==tempStatus1Muon->pdgId() ){
tempStatus1Muon = tempStatus1Muon->daughter(i);
status = tempStatus1Muon->status();
for (unsigned int i=0; i<tempMuon->numberOfDaughters(); ++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<const reco::Candidate*> daughters;
for (unsigned int i=0; i<tempMuon->numberOfDaughters(); ++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)
{
Expand Down Expand Up @@ -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;
}
}
}
Expand All @@ -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 = "<<motherPdgId<<std::endl;
Expand All @@ -430,6 +458,7 @@ GenMuonPair MuScleFitMuonSelector::findGenMuFromRes( const reco::GenParticleColl
// muFromRes.second = (lorentzVector(status1Muon->p4().px(),status1Muon->p4().py(),
// status1Muon->p4().pz(),status1Muon->p4().e()));
}
muFromRes.motherId = motherPdgId;
}
}
}
Expand All @@ -455,7 +484,7 @@ std::pair<lorentzVector, lorentzVector> 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)
Expand Down
Expand Up @@ -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<MuScleFitMuon> & muons,
Expand Down

0 comments on commit 1430449

Please sign in to comment.