Merge pull request #10477 from emiglior/changes_pythia8_74X

Update for GEN infos in 74X or later

MuonAnalysis/MomentumScaleCalibration/bin/TreeFromDump.cc

@@ -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];

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

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

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)
 

MuonAnalysis/MomentumScaleCalibration/plugins/BuildFile.xml

@@ -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"/>

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<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 ) {
@@ -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!
 				  );
 
 

MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitMuonSelector.cc

@@ -8,24 +8,45 @@ 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)
 {
   reco::TrackRef iTrack = aMuon->innerTrack();
@@ -217,7 +238,9 @@ void MuScleFitMuonSelector::selectMuons(const edm::Event & event, std::vector<Mu
   if( speedup_ == false )
   {
     if( PATmuons_ ) {
-      selectGeneratedMuons(collAll, collMuSel, genPair, plotter);
+      // EM 2015.04.02 temporary fix to run on MINIAODSIM (which contains PAT muons) but not the "onia2MuMuPatTrkTrk" collection   
+      // selectGeneratedMuons(collAll, collMuSel, genPair, plotter);
+      selectGenSimMuons(event, genPair, simPair, plotter);
     }
     else {
       selectGenSimMuons(event, genPair, simPair, plotter);
@@ -384,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; 
       }
     }
   }
@@ -408,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;
@@ -428,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; 
       }
     }
   }
@@ -453,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)

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<MuScleFitMuon> & muons,