Skip to content

Instantly share code, notes, and snippets.

@mmusich

mmusich/MakeFinalMuoTreeMINIAOD.cc

Created February 4, 2024 16:56
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mmusich/fd03bc1384c07f69fce58c835aaec568 to your computer and use it in GitHub Desktop.
Save mmusich/fd03bc1384c07f69fce58c835aaec568 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
MakeFinalMuoTreeMINIAOD.cc
#include <memory>
#include <iostream>
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "FWCore/Common/interface/TriggerResultsByName.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/Candidate/interface/VertexCompositeCandidate.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/Common/interface/ValueMap.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
//#include "DataFormats/PatCandidates/interface/Electron.h"
//#include "DataFormats/PatCandidates/interface/VIDCutFlowResult.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/MET.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
//#include "DataFormats/PatCandidates/interface/IsolatedTrack.h"
#include "DataFormats/PatCandidates/interface/PackedTriggerPrescales.h"
#include "DataFormats/PatCandidates/interface/PackedGenParticle.h"
#include "DataFormats/METReco/interface/PFMET.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/TrackReco/interface/DeDxHitInfo.h"
#include "DataFormats/TrackReco/interface/DeDxData.h"
#include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
//#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
//#include "DataFormats/EgammaCandidates/interface/Photon.h"
//#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
#include "DataFormats/GeometrySurface/interface/BoundDisk.h"
#include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/Records/interface/CaloTopologyRecord.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
//#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
//#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
#include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
#include "TrackingTools/TrackAssociator/interface/FiducialVolume.h"
#include "TrackingTools/TrackAssociator/interface/DetIdInfo.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
//#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
//#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "RecoTracker/DeDx/interface/DeDxTools.h"
#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
//#include "TH1F.h"
//#include "TH1S.h"
#include "TClonesArray.h"
#include "TObject.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TLorentzVector.h"
#include "TMath.h"
using namespace std;
typedef math::XYZTLorentzVector LorentzVector;
typedef math::XYZPoint Point;
class MakeFinalMuoTreeMINIAOD : public edm::EDAnalyzer {
public:
MakeFinalMuoTreeMINIAOD(const edm::ParameterSet&);
~MakeFinalMuoTreeMINIAOD() {};
static void fillDescriptions(edm::ConfigurationDescriptions&);
private:
virtual void analyze(const edm::Event&, const edm::EventSetup&);
virtual void endJob();
virtual void beginRun(const edm::Run& , const edm::EventSetup&);
TLorentzVector lorentzMomentum(const LorentzVector) const;
TLorentzVector lorentzMomentum(double, double, double, double) const;
bool refitPV(const edm::Event&, const edm::EventSetup&);
static const unsigned int maxAccTrig=500;
static const unsigned int maxLep=2000;
//const double eleMass = 0.000510998950;
//const double piMass = 0.13957018;
//const double protonMass = 0.938272046;
//const double lambdaMass = 1.115683;
const double muMass = 0.1056583745;
const int minHits = 5;
const float MeVperADCStrip = 3.61e-06*265;
const float cutMinCL = 0;
std::string theOutFileName;
std::string era;
TFile *outFile;
TTree *fTree;
unsigned int eventNb;
unsigned int runNb;
unsigned int lumiBlock;
int nvtx;
float nputrue;
float genweight;
float genDiLepMass;
int genLepCode;
bool metFilters;
double prefiringweightECAL_;
double prefiringweightECALup_;
double prefiringweightECALdown_;
double prefiringweightMUON_;
double prefiringweightMUONup_;
double prefiringweightMUONdown_;
float pfMETx;
float pfMETy;
float pfMET;
float ht;
int iAcceptedTriggers_;
unsigned int nGen_;
int idGen_[maxLep];
unsigned int nMuo_;
int idMuo_[maxLep];
int nlMuo_[maxLep];
int fHLTMuo_[maxLep];
float isoMuoTR_[maxLep];
float isoMuoPF_[maxLep];
float sip3dMuo_[maxLep];
float dxyMuo_[maxLep];
float dzMuo_[maxLep];
float dedxMuo_[maxLep];
float ptErrorMuo_[maxLep];
int dedxMuoNhits_[maxLep];
double clVertexMuo2_;
double L_sigmaMuo2_;
double L_absMuo2_;
double L_3dDistMuo2_;
double cosphiMuo2_;
double cosphiMuoOri2_;
double muo0dis_, muo1dis_;
int codeVertexMuo0_;
int codeVertexMuo1_;
double clVertexMuo2_sys_;
double L_sigmaMuo2_sys_;
double L_absMuo2_sys_;
double L_3dDistMuo2_sys_;
double cosphiMuo2_sys_;
double deltaSecVertPrimVertMuo2_[3];
double PrimVertGen_[3];
double PrimVertRec_[3];
double PrimVertRefRec_[3];
double SecVertGen_[3];
double SecVertMuoRec_[3];
float dimuoMass_;
int dimuoCharge_;
double deltaSecVertPrimVertMuo2_sys_[3];
double SecVertMuoRec_sys_[3];
float dimuoMass_sys_;
// for muons dedx
Float_t muodedxArray[100];
Short_t muoshapeArray[100];
TClonesArray* gen;
TClonesArray* muo;
TClonesArray* muotrack;
TClonesArray* acceptedTriggers;
const edm::EDGetTokenT<GenEventInfoProduct> genEvtInfoToken_;
const edm::EDGetTokenT<std::vector<pat::PackedGenParticle> > genParticleToken_;
const edm::EDGetTokenT<edm::View<reco::GenParticle> > mcTruthCollectionToken_; // prunedGenParticles
const edm::EDGetTokenT<std::vector<PileupSummaryInfo> > puSummaryToken_;
const edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
const edm::EDGetTokenT<std::vector<pat::MET> > metToken_;
const edm::EDGetTokenT<std::vector<pat::Muon> > muonToken_;
const edm::EDGetTokenT<edm::TriggerResults> metFiltersBitToken_;
const edm::EDGetTokenT<edm::TriggerResults> triggerBitToken_;
const edm::EDGetTokenT<std::vector<pat::PackedCandidate> > pfCandToken_;
const edm::EDGetTokenT<reco::DeDxHitInfoAss> dedxMapToken_;
const edm::EDGetTokenT<std::vector<reco::DeDxHitInfo> > dedxHitInfoToken_;
const edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackRecordToken_;
edm::EDGetTokenT< double > prefweight_token;
edm::EDGetTokenT< double > prefweightup_token;
edm::EDGetTokenT< double > prefweightdown_token;
edm::EDGetTokenT< double > prefweightECAL_token;
edm::EDGetTokenT< double > prefweightECALup_token;
edm::EDGetTokenT< double > prefweightECALdown_token;
edm::EDGetTokenT< double > prefweightMUON_token;
edm::EDGetTokenT< double > prefweightMUONup_token;
edm::EDGetTokenT< double > prefweightMUONdown_token;
edm::EDGetTokenT< bool > BadPFMuonDzFilterToken_;
const edm::TriggerNames * metFiltersNames_;
const edm::TriggerNames * triggerNames_;
HLTPrescaleProvider hltPrescaleProvider_;
// get stuff from Event Setup
//edm::ESHandle<TransientTrackBuilder> theB;
public:
int sysRun_;
edm::Handle<GenEventInfoProduct> genEvtInfo_;
edm::Handle<std::vector<pat::PackedGenParticle> > genParticles_;
edm::Handle<edm::View<reco::GenParticle> > mcTruthCollection_;
edm::Handle<std::vector<PileupSummaryInfo> > puSummaryInfos_;
edm::Handle<reco::VertexCollection> vertices_;
edm::Handle<std::vector<pat::MET> > mets_;
edm::Handle<std::vector<pat::Muon> > muons_;
edm::Handle<edm::TriggerResults> metFilters_;
edm::Handle<edm::TriggerResults> triggerBits_;
edm::Handle<std::vector<pat::PackedCandidate> > pfcands_;
edm::Handle<reco::DeDxHitInfoAss> dedxMap_;
edm::Handle<std::vector<reco::DeDxHitInfo> > dedxHitInfo_;
edm::Handle<reco::BeamSpot> beamSpot_;
edm::Handle<bool> BadPFMuonDzFilterHandle_;
edm::Handle< double > theprefweight;
edm::Handle< double > theprefweightup;
edm::Handle< double > theprefweightdown;
edm::Handle< double > theprefweightECAL;
edm::Handle< double > theprefweightECALup;
edm::Handle< double > theprefweightECALdown;
edm::Handle< double > theprefweightMUON;
edm::Handle< double > theprefweightMUONup;
edm::Handle< double > theprefweightMUONdown;
};
void MakeFinalMuoTreeMINIAOD::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
eventNb = iEvent.id().event();
runNb = iEvent.id().run();
lumiBlock = iEvent.luminosityBlock();
//cout << "*************************************************" << eventNb << " " << runNb << " " << lumiBlock << "*****************" << endl;
vector<reco::TransientTrack> tks;
vector<reco::TransientTrack> tks_sys;
iEvent.getByToken(vtxToken_, vertices_);
iEvent.getByToken(metFiltersBitToken_, metFilters_);
iEvent.getByToken(triggerBitToken_, triggerBits_);
iEvent.getByToken(muonToken_,muons_);
iEvent.getByToken(metToken_, mets_);
iEvent.getByToken(pfCandToken_, pfcands_);
iEvent.getByToken(dedxMapToken_,dedxMap_);
iEvent.getByToken(dedxHitInfoToken_,dedxHitInfo_);
iEvent.getByToken(beamSpotToken_, beamSpot_);
iEvent.getByToken(BadPFMuonDzFilterToken_,BadPFMuonDzFilterHandle_);
// get the transient track builder
auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
// GEN Information
genweight = 1;
genLepCode = 0;
genDiLepMass = 0;
nGen_ = 0;
prefiringweightECAL_ = 1;
prefiringweightECALup_ = 1;
prefiringweightECALdown_ = 1;
prefiringweightMUON_ = 1;
prefiringweightMUONup_ = 1;
prefiringweightMUONdown_ = 1;
PrimVertGen_[0] =0;PrimVertGen_[1] =0;PrimVertGen_[2] =0;
PrimVertRec_[0] =0;PrimVertRec_[1] =0;PrimVertRec_[2] =0;
SecVertGen_[0]=0;SecVertGen_[1]=0;SecVertGen_[2]=0;
SecVertMuoRec_[0]=0;SecVertMuoRec_[1]=0;SecVertMuoRec_[2]=0;
if( !iEvent.isRealData() ) {
iEvent.getByToken(genEvtInfoToken_, genEvtInfo_);
iEvent.getByToken(genParticleToken_, genParticles_);
iEvent.getByToken(mcTruthCollectionToken_, mcTruthCollection_);
iEvent.getByToken(puSummaryToken_, puSummaryInfos_);
genweight = genEvtInfo_->weight();
//cout << "Compute preFire weight" << endl;
iEvent.getByToken(prefweight_token, theprefweight ) ;
iEvent.getByToken(prefweightup_token, theprefweightup ) ;
iEvent.getByToken(prefweightdown_token, theprefweightdown ) ;
iEvent.getByToken(prefweightECAL_token, theprefweightECAL ) ;
iEvent.getByToken(prefweightECALup_token, theprefweightECALup ) ;
iEvent.getByToken(prefweightECALdown_token, theprefweightECALdown ) ;
iEvent.getByToken(prefweightMUON_token, theprefweightMUON ) ;
iEvent.getByToken(prefweightMUONup_token, theprefweightMUONup ) ;
iEvent.getByToken(prefweightMUONdown_token, theprefweightMUONdown ) ;
//cout << "Prefiring: " << (*theprefweight) << " " << (*theprefweightECAL)*(*theprefweightMUON) << " | " << (*theprefweightECAL) << " " << (*theprefweightMUON) << endl;
prefiringweightECALup_ = (*theprefweightECALup);
prefiringweightECALdown_ = (*theprefweightECALdown);
prefiringweightECAL_ = (*theprefweightECAL);
prefiringweightECALup_ = (*theprefweightECALup);
prefiringweightECALdown_ = (*theprefweightECALdown);
prefiringweightMUON_ = (*theprefweightMUON);
prefiringweightMUONup_ = (*theprefweightMUONup);
prefiringweightMUONdown_ = (*theprefweightMUONdown);
Point theVertex;
int idPart = 0;
bool first = true;
bool goodMCvertex = false;
for (size_t i = 0; i < mcTruthCollection_->size(); i++) {
if( abs((*mcTruthCollection_)[i].pdgId()) == 11 ||
abs((*mcTruthCollection_)[i].pdgId()) == 13 ) {
if ( (*mcTruthCollection_)[i].isPromptFinalState() && first) {
//cout << "first pruned prompt " << i << " " << (*mcTruthCollection_)[i].pdgId() << " " << (*mcTruthCollection_)[i].vertex() << " " << (*mcTruthCollection_)[i].mother() << endl;
theVertex = (*mcTruthCollection_)[i].vertex();
idPart = (*mcTruthCollection_)[i].pdgId();
first = false;
}
if ( (*mcTruthCollection_)[i].isPromptFinalState() && !first &&
(*mcTruthCollection_)[i].vertex() == theVertex &&
abs( (*mcTruthCollection_)[i].pdgId() ) == abs(idPart) &&
(*mcTruthCollection_)[i].pdgId()*idPart < 0 ) {
//cout << "second pruned prompt " << i << " " << (*mcTruthCollection_)[i].pdgId() << " " << (*mcTruthCollection_)[i].vertex() << " " << (*mcTruthCollection_)[i].mother() << endl;
goodMCvertex = true;
}
//if ( (*mcTruthCollection_)[i].isDirectPromptTauDecayProductFinalState() )
// cout << "pruned fromtau " << i << " " << (*mcTruthCollection_)[i].pdgId() << " " << (*mcTruthCollection_)[i].vertex() << " " << (*mcTruthCollection_)[i].mother() << endl;
}
if( (*mcTruthCollection_)[i].pdgId()==23 || (*mcTruthCollection_)[i].pdgId()==32) {
//cout << "Z0/P: " << i << " " << (*mcTruthCollection_)[i].pdgId() << " " << (*mcTruthCollection_)[i].vertex() << " " << (*mcTruthCollection_)[i].mother() << endl;
for (unsigned idx = 0; idx < (*mcTruthCollection_)[i].numberOfDaughters(); ++idx) {
//cout << " dau: " << (*mcTruthCollection_)[i].daughter(idx)->pdgId() << endl;
int idDau = abs((*mcTruthCollection_)[i].daughter(idx)->pdgId());
if( idDau==11 || idDau==13) {
//cout << "Prim Vertex Gen: " << (*mcTruthCollection_)[i].vertex() << endl;
PrimVertGen_[0]=(*mcTruthCollection_)[i].vx();
PrimVertGen_[1]=(*mcTruthCollection_)[i].vy();
PrimVertGen_[2]=(*mcTruthCollection_)[i].vz();
}
}
if( (*mcTruthCollection_)[i].mother(0)->pdgId() !=23 && (*mcTruthCollection_)[i].mother(0)->pdgId() !=32 ) {
PrimVertGen_[0]=(*mcTruthCollection_)[i].vx();
PrimVertGen_[1]=(*mcTruthCollection_)[i].vy();
PrimVertGen_[2]=(*mcTruthCollection_)[i].vz();
//cout << "pv " << PrimVertGen_[0] << " " << PrimVertGen_[1] << " " << PrimVertGen_[2] << endl;
}
if( (*mcTruthCollection_)[i].daughter(0)->pdgId() !=23 && (*mcTruthCollection_)[i].daughter(0)->pdgId() !=32 ) {
SecVertGen_[0]=(*mcTruthCollection_)[i].daughter(0)->vx();
SecVertGen_[1]=(*mcTruthCollection_)[i].daughter(0)->vy();
SecVertGen_[2]=(*mcTruthCollection_)[i].daughter(0)->vz();
//cout << "sv " << SecVertGen_[0] << " " << SecVertGen_[1] << " " << SecVertGen_[2] << endl;
}
}
}
if (goodMCvertex) {
//cout << "********* GOOD GEN VERTEX *********" << theVertex << endl;
//cout << "Pri Vertex Gen: " << PrimVertGen_[0] << " " << PrimVertGen_[1] << " " << PrimVertGen_[2] << endl;
//cout << "Sec Vertex Gen: " << theVertex << endl;
SecVertGen_[0]=theVertex.x();
SecVertGen_[1]=theVertex.y();
SecVertGen_[2]=theVertex.z();
}
int nEle=0;
int nMuo=0;
//int nTau=0;
//cout << "evt wgt:" << genEvtInfo_->weight() << endl;
//TLorentzVector sum( 0, 0, 0, 0 );
for(unsigned int iP = 0; iP < genParticles_->size(); ++iP){
const pat::PackedGenParticle * part = &genParticles_->at(iP);
//if( abs(part->pdgId())==15 ) cout << "allParticle: " << iP << " " << part->pdgId() << " pt/eta/phi: " << part->pt() << " " << part->eta() << " " << part->phi() <<
// " px/py/pz/E: " << part->px() << " " << part->py() << " " << part->pz() << " " << part->energy() <<
// " vertex: " << part->vx() << " " << part->vy() << " " << part->vz() << endl;
//reco::GenParticleRef genRef = reco::GenParticleRef(genParticles_, iP);
if( part->isPromptFinalState() ) {
//if ( abs(part->pdgId())==19 ) cout << "NO: " << part->pdgId() << " " << part->status() << " " << part->mass() << " " << part->charge() << " "
// << part->pt() << " " << part->eta() << " " << part->phi() << endl;
//TLorentzVector p4( lorentzMomentum( part->p4() ) );
if( abs(part->pdgId() ) == 11 ) {
//cout << "ele found: " << nMuo << " pt/eta/phi: " << part->pt() << " " << part->eta() << " " << part->phi() <<
// " px/py/pz/E: " << part->px() << " " << part->py() << " " << part->pz() << " " << part->energy() <<
// " vertex: " << part->vx() << " " << part->vy() << " " << part->vz() << endl;
idGen_[nGen_] = part->pdgId();
new((*gen)[nGen_]) TLorentzVector( lorentzMomentum( part->p4() ) );
nEle++;
nGen_++;
//sum += p4;
}
if( abs(part->pdgId() ) == 13 ) {
//cout << "muon found: " << nMuo << " pt/eta/phi: " << part->pt() << " " << part->eta() << " " << part->phi() <<
// " px/py/pz/E: " << part->px() << " " << part->py() << " " << part->pz() << " " << part->energy() <<
// " vertex: " << part->vx() << " " << part->vy() << " " << part->vz() << endl;
idGen_[nGen_] = part->pdgId();
new((*gen)[nGen_]) TLorentzVector( lorentzMomentum( part->p4() ) );
if( part->isDirectPromptTauDecayProductFinalState() ) cout << "prompt muo from tau" << endl;
if( part->isDirectHardProcessTauDecayProductFinalState() ) cout << "hard muo from tau" << endl;
nMuo++;
nGen_++;
//sum += p4;
}
}
if( ( part->isDirectPromptTauDecayProductFinalState() || part->isDirectHardProcessTauDecayProductFinalState() ) &&
abs(part->pdgId() ) == 11 ) {
idGen_[nGen_] = 1000011*part->pdgId()/abs(part->pdgId());
new((*gen)[nGen_]) TLorentzVector( lorentzMomentum( part->p4() ) );
nEle++;
nGen_++;
}
if( ( part->isDirectPromptTauDecayProductFinalState() || part->isDirectHardProcessTauDecayProductFinalState() ) &&
abs(part->pdgId() ) == 13 ) {
idGen_[nGen_] = 1000013*part->pdgId()/abs(part->pdgId());
new((*gen)[nGen_]) TLorentzVector( lorentzMomentum( part->p4() ) );
nMuo++;
nGen_++;
}
}
//cout << "GEN Particles: nEle: " << nEle << " nMuo: " << nMuo << endl;
}
triggerNames_ = &iEvent.triggerNames(*triggerBits_);
metFiltersNames_ = &iEvent.triggerNames(*metFilters_);
// Pileup Reweighting
nputrue = 0;
if(puSummaryInfos_.isValid()) {
for( const auto& psu : *puSummaryInfos_ ) {
if(psu.getBunchCrossing() == 0) nputrue = psu.getTrueNumInteractions();
}
}
/* // CHECK FOR ZPrime VERTEX
float minDist = 1000;
int idxNearVert = -1;
for(unsigned int nVtx=0; nVtx < vertices_->size(); nVtx++) {
reco::Vertex pVtx = (*vertices_)[nVtx];
float distGenPV = sqrt(pow(PrimVertGen_[0]-pVtx.position().x(),2)+pow(PrimVertGen_[1]-pVtx.position().y(),2)+pow(PrimVertGen_[2]-pVtx.position().z(),2));
if (distGenPV<minDist) {
minDist=distGenPV;
idxNearVert=nVtx;
//cout << "PV FOUND: " << nVtx << " | " << pVtx.position() << " | " << distGenPV << endl;
}
}
if(idxNearVert>=0) cout << "PV FOUND: " << idxNearVert << " | " << (*vertices_)[idxNearVert].position() << " | " << minDist << endl;
*/
// Vertex INFO
nvtx = vertices_->size();
int primaryVertexIndex_;
double cov_pv[3][3];
if(vertices_->size() > 0) primaryVertexIndex_ = 0;
//bool hasgoodvtx = false;
//if(!(*vertices_)[primaryVertexIndex_].isFake() &&
// (*vertices_)[primaryVertexIndex_].ndof() > 4. &&
// (*vertices_)[primaryVertexIndex_].position().Rho() <= 2.0 &&
// fabs((*vertices_)[primaryVertexIndex_].position().Z())<=24.0) hasgoodvtx = true;
reco::Vertex primaryVertex = (*vertices_)[primaryVertexIndex_];
//cout << "Prim Vertex RECO: " << primaryVertex.position() << endl;
PrimVertRec_[0]=primaryVertex.position().x();
PrimVertRec_[1]=primaryVertex.position().y();
PrimVertRec_[2]=primaryVertex.position().z();
cov_pv[0][0] = primaryVertex.covariance(0,0);
cov_pv[1][0] = primaryVertex.covariance(1,0);
cov_pv[2][0] = primaryVertex.covariance(2,0);
cov_pv[0][1] = cov_pv[1][0];
cov_pv[1][1] = primaryVertex.covariance(1,1);
cov_pv[2][1] = primaryVertex.covariance(2,1);
cov_pv[0][2] = cov_pv[2][0];
cov_pv[1][2] = cov_pv[2][1];
cov_pv[2][2] = primaryVertex.covariance(2,2);
//MET Filters
// https://twiki.cern.ch/twiki/bin/viewauth/CMS/MissingETOptionalFiltersRun2#UL_data
metFilters = true;
for(unsigned int i = 0; i < metFilters_->size(); ++i) {
const auto &trigname = metFiltersNames_->triggerName(i);
if(trigname=="Flag_goodVertices" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_globalSuperTightHalo2016Filter" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_EcalDeadCellTriggerPrimitiveFilter" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_BadPFMuonFilter" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_BadPFMuonDzFilter" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_eeBadScFilter" && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_hfNoisyHitsFilter" && (era=="2022" ) && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_HBHENoiseFilter" && (era=="2016" || era=="2017" || era=="2018") && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_HBHENoiseIsoFilter" && (era=="2016" || era=="2017" || era=="2018") && !metFilters_->accept(i) ) metFilters=false;
if(trigname=="Flag_ecalBadCalibFilter" && (era=="2017" || era=="2018" || era=="2022") && !metFilters_->accept(i) ) metFilters=false;
}
// FILTRO PRESENTE SOLO IN MiniAODv2 -> Se MiniAODv1 lo calcola qui
//cout << "metFiltersB: " << metFilters << endl;
//cout << "BadPFMuonDzFilterHandle_ " << *BadPFMuonDzFilterHandle_ << endl;
if (!(*BadPFMuonDzFilterHandle_)) metFilters=false;
//cout << "metFiltersA: " << metFilters << endl;
// FIRED TRIGGERS
HLTConfigProvider const& hltConfig = hltPrescaleProvider_.hltConfigProvider();
std::vector<std::string> activeHLTPathsInThisEvent = hltConfig.triggerNames();
iAcceptedTriggers_ = 0 ;
for (std::vector<std::string>::const_iterator iHLT = activeHLTPathsInThisEvent.begin();
iHLT != activeHLTPathsInThisEvent.end();
++iHLT) {
//matching with regexp filter name. More than 1 matching filter is allowed
TObjString *str= (TObjString*)acceptedTriggers->ConstructedAt(iAcceptedTriggers_);
str->SetString(iHLT->c_str());
//auto const prescales = hltPrescaleProvider_.prescaleValues(iEvent, iSetup, *iHLT);
//cout << prescales.first << " " << prescales.second << " " << *iHLT << endl;
++iAcceptedTriggers_;
}
/*
iAcceptedTriggers_ = 0 ;
for(unsigned int i = 0; i < triggerBits_->size(); ++i) {
const auto &trigname = triggerNames_->triggerName(i);
if( triggerBits_->accept(i) ) {
const std::pair<int,int> prescales(hltConfig_.prescaleValues(iEvent,iSetup,trigname));
TObjString *str= (TObjString*)acceptedTriggers->ConstructedAt(iAcceptedTriggers_);
str->SetString(trigname.data());
//cout << trigname.data() << endl;
++iAcceptedTriggers_;
}
}
*/
// get PFMet
pfMET=-1;
pfMETx=99999;
pfMETy=99999;
const std::vector<pat::MET> *metCol = mets_.product();
const pat::MET theMet = metCol->front();
pfMET = theMet.et();
pfMETx = theMet.px();
pfMETy = theMet.py();
//=========================//
// MUONS //
//=========================//
// Variables for Kalman filter
clVertexMuo2_ = -1;
clVertexMuo2_sys_ = -1;
L_absMuo2_ = -1;
L_absMuo2_sys_ = -1;
L_3dDistMuo2_ = 9999;
L_3dDistMuo2_sys_ = 9999;
L_sigmaMuo2_ = 0;
L_sigmaMuo2_sys_ = 0;
cosphiMuo2_ = 9999;
cosphiMuo2_sys_ = 9999;
cosphiMuoOri2_ = 9999;
codeVertexMuo0_ = -1;
codeVertexMuo1_ = -1;
muo0dis_ = -1;
muo1dis_ = -1;
dimuoMass_ = 0;
dimuoMass_sys_ = 0;
dimuoCharge_ = 0;
double derivMuo[3];
double derivMuo_sys[3];
double cov_svMuo2[3][3];
double cov_svMuo2_sys[3][3];
double pxDiMuo = 0, pyDiMuo = 0, pzDiMuo = 0, piDiMuo = 0,
pxDiMuoOri = 0, pyDiMuoOri = 0, pzDiMuoOri = 0, piDiMuoOri = 0;
double pxDiMuo_sys = 0, pyDiMuo_sys = 0, pzDiMuo_sys = 0, piDiMuo_sys = 0;
//double dimuoMassOri_ = 0;
deltaSecVertPrimVertMuo2_[0]=0;deltaSecVertPrimVertMuo2_[1]=0;deltaSecVertPrimVertMuo2_[2]=0;
deltaSecVertPrimVertMuo2_sys_[0]=0;deltaSecVertPrimVertMuo2_sys_[1]=0;deltaSecVertPrimVertMuo2_sys_[2]=0;
TransientVertex aTransientVertexMuo01;
TransientVertex aTransientVertexMuo01_sys;
GlobalPoint theSecondaryVertexPositionMuo;
GlobalPoint theSecondaryVertexPositionMuo_sys;
const int nTotMuonHLTpath = 15;
char const *hltstringMuo[nTotMuonHLTpath] = {
"HLT_IsoMu24_v*",
"HLT_IsoMu27_v*",
"HLT_IsoTkMu24_v*",
"HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_v*",
"HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_v*",
"HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_v*",
"HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ_v*",
"HLT_TkMu17_TrkIsoVVL_TkMu8_TrkIsoVVL_v*",
"HLT_TkMu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ_v*",
"HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass3p8_v*",
"HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass8_v*",
"HLT_Mu50_v*",
"HLT_TkMu50_v*",
"HLT_OldMu100_v*",
"HLT_TkMu100_v*"
};
//vector<pat::Muon> myGoodMuonVector;
vector<const reco::Candidate *> myGoodMuonVector;
vector<const reco::Track *> myGoodTrackMuonVector;
vector<int> goodMuons;
reco::TransientTrack muon0, muon1;
reco::TransientTrack muon0_sys, muon1_sys;
reco::TransientTrack muon0refitted, muon1refitted;
reco::TransientTrack muon0refitted_sys, muon1refitted_sys;
TransientVertex theSecondaryVertexMuo, theSecondaryVertexMuo_sys;
nMuo_ = 0;
for (const pat::Muon &muon : *muons_) {
//if (muon.pt()<5) continue;
idMuo_[nMuo_] = 1;
nlMuo_[nMuo_] = -1;
fHLTMuo_[nMuo_] = 0;
dxyMuo_[nMuo_] = -1;
dzMuo_[nMuo_] = 9999;
isoMuoPF_[nMuo_] = -1;
isoMuoTR_[nMuo_] = -1;
dedxMuo_[nMuo_] = -1;
dedxMuoNhits_[nMuo_] = -1;
ptErrorMuo_[nMuo_] = 9999;
//cout << muon.muonBestTrack().get() << " " << muon.innerTrack().get() << " " << muon.combinedMuon().get() << " " << muon.globalTrack().get() << endl;
const reco::TrackRef t = muon.innerTrack();
//const reco::TrackRef t = muon.muonBestTrack();
if ( ! t.isNull() ) {
nlMuo_[nMuo_] = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement();
dxyMuo_[nMuo_] = muon.innerTrack()->dxy(primaryVertex.position());
dzMuo_[nMuo_] = muon.innerTrack()->dz(primaryVertex.position());
//nlMuo_[nMuo_] = muon.muonBestTrack()->hitPattern().trackerLayersWithMeasurement();
//dxyMuo_[nMuo_] = muon.muonBestTrack()->dxy(primaryVertex.position());
//dzMuo_[nMuo_] = muon.muonBestTrack()->dz(primaryVertex.position());
if(t->quality(reco::TrackBase::highPurity) ) {
//cout << "innerTrack: " << t->pt() << " " << t->eta() << " " << t->phi() << endl;
if(
t->chi2()/t->ndof() <= 2.5 &&
t->numberOfValidHits() >= 5 &&
t->hitPattern().numberOfValidPixelHits() >= 2
//fabs(t->dxy(primaryVertex.position())) < 0.1 &&
//fabs(t->dz(primaryVertex.position())) < 1
) {
myGoodMuonVector.push_back(&muon);
myGoodTrackMuonVector.push_back( &*t );
goodMuons.push_back(nMuo_);
}
}
}
//cout << "myGoodTrackMuonVector : " << myGoodTrackMuonVector.size() << endl;
for(int nHLT=0;nHLT<nTotMuonHLTpath;nHLT++) {
if( muon.triggered( hltstringMuo[nHLT] ) ) {
fHLTMuo_[nMuo_] |= 1UL << nHLT;
//cout << "** " << nHLT << " " << hltstringMuo[nHLT] << " " << fHLTMuo_[nMuo_] << endl;
}
}
//cout << "----- " << nMuo_ << " ftr " << fHLTMuo_[nMuo_] << " " << muon.pt() << endl;
// SEE https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideMuonIdRun2#Muon_Isolation
if (muon.isLooseMuon()) idMuo_[nMuo_] +=10 ;
if (muon.isMediumMuon()) idMuo_[nMuo_] +=100 ;
if( muon.isGlobalMuon() &&
muon.isPFMuon() &&
muon.globalTrack()->normalizedChi2() < 10. &&
muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0 &&
muon.numberOfMatchedStations() > 1 &&
muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0 &&
muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 ) idMuo_[nMuo_] +=1000 ;
if (muon.isTightMuon(primaryVertex)) idMuo_[nMuo_] +=10000 ;
if (muon.isSoftMuon(primaryVertex)) idMuo_[nMuo_] +=100000 ;
if (muon.isHighPtMuon(primaryVertex)) idMuo_[nMuo_] +=1000000 ;
if (muon.passed( reco::Muon::InTimeMuon )) idMuo_[nMuo_] +=10000000 ;
isoMuoPF_[nMuo_] = ( muon.pfIsolationR04().sumChargedHadronPt + max(0.0 , muon.pfIsolationR04().sumNeutralHadronEt + muon.pfIsolationR04().sumPhotonEt - 0.5 * muon.pfIsolationR04().sumPUPt) ) /muon.pt();
isoMuoTR_[nMuo_] = muon.isolationR03().sumPt/muon.pt();
//cout << "ISO " << (isoMuoPF_[nMuo_]<0.15) << " " << muon.passed(reco::Muon::PFIsoTight) << " " << isoMuoPF_[nMuo_] << endl;
//cout << "XTC " << muon.pt() << " " << muon.p() << " " << muon.energy() << " " << muon.eta() << " " << muon.phi() << " " << muon.mass() << endl;
new((*muo)[nMuo_]) TLorentzVector( lorentzMomentum( muon.p4() ) );
idMuo_[nMuo_] *= muon.charge();
float sip3d=fabs(muon.dB(muon.PV3D) / muon.edB(muon.PV3D)); // https://github.com/cms-sw/cmssw/blob/2c1d10a3c57a4be2640000233fcf5fb8b8b158bb/PhysicsTools/NanoAOD/python/electrons_cff.py#L373
sip3dMuo_[nMuo_] = sip3d;
if (muon.isGlobalMuon() ) ptErrorMuo_[nMuo_] = muon.globalTrack()->ptError();
else if (muon.isTrackerMuon()) ptErrorMuo_[nMuo_] = muon.innerTrack()->ptError();
nMuo_++;
}
//if(goodMuons.size()>1 && goodMuons[0]==0 && goodMuons[1]==1) cout << "--OUT "<< fHLTMuo_[0] << " END--" << endl;
//for(unsigned int n=0;n<myGoodMuonVector.size();n++) cout << "mgvm: " << myGoodMuonVector[n]->pt() << " " << myGoodMuonVector[n]->charge() << endl;
//cout << "CHECK NMUONS: " << nMuoNOref_ << " " << nMuoSIref_++ << endl;
if( myGoodMuonVector.size() >= 2 ) {
dimuoCharge_ = myGoodMuonVector[0]->charge() * myGoodMuonVector[1]->charge();
const reco::Track *t0 = myGoodTrackMuonVector[0];
const reco::Track *t1 = myGoodTrackMuonVector[1];
new((*muotrack)[0]) TLorentzVector( t0->px(), t0->py(), t0->pz(), t0->p() );
new((*muotrack)[1]) TLorentzVector( t1->px(), t1->py(), t1->pz(), t1->p() );
double d00new = t0->d0();
double dz0new = t0->dz();
double phi0new = t0->phi();
double theta0new = t0->theta();
double pt0new = t0->pt();
double d01new = t1->d0();
double dz1new = t1->dz();
double phi1new = t1->phi();
double theta1new = t1->theta();
double pt1new = t1->pt();
switch (sysRun_) {
case 0:
// NO SYS
break;
case 1:
d00new = t0->d0() + t0->d0Error();
d01new = t1->d0() + t1->d0Error();
break;
case -1:
d00new = t0->d0() - t0->d0Error();
d01new = t1->d0() - t1->d0Error();
break;
case 2:
dz0new = t0->dz() + t0->dzError();
dz1new = t1->dz() + t1->dzError();
break;
case -2:
dz0new = t0->dz() - t0->dzError();
dz1new = t1->dz() - t1->dzError();
break;
case 3:
theta0new = t0->theta() + t0->thetaError();
theta1new = t1->theta() + t1->thetaError();
break;
case -3:
theta0new = t0->theta() - t0->thetaError();
theta1new = t1->theta() - t1->thetaError();
break;
case 4:
phi0new = t0->phi() + t0->phiError();
phi1new = t1->phi() + t1->phiError();
break;
case -4:
phi0new = t0->phi() - t0->phiError();
phi1new = t1->phi() - t1->phiError();
break;
case 5:
pt0new = t0->pt() + t0->ptError();
pt1new = t1->pt() + t1->ptError();
break;
case -5:
pt0new = t0->pt() - t0->ptError();
pt1new = t1->pt() - t1->ptError();
break;
default:
cout << "ERROR IN sysRun_ : " << sysRun_ << endl;
}
double px0new = pt0new*cos(phi0new);
double py0new = pt0new*sin(phi0new);
double pz0new = pt0new/tan(theta0new);
double vz0new = t0->referencePoint().Z();
double vy0new = (vz0new-dz0new)*py0new/pz0new - d00new*px0new/pt0new;
double vx0new = d00new*pt0new/py0new + vy0new*px0new/py0new;
math::XYZPoint refPoint0NEW(vx0new,vy0new,vz0new);
math::XYZVector momentum0NEW(px0new,py0new,pz0new);
const reco::Track *t0_sys = new reco::Track(
t0->chi2(),
t0->ndof(),
refPoint0NEW,
momentum0NEW,
t0->charge(),
t0->covariance(),
t0->algo(),
reco::TrackBase::highPurity,
t0->t0(), t0->beta(), t0->covt0t0(), t0->covBetaBeta()
);
double px1new = pt1new*cos(phi1new);
double py1new = pt1new*sin(phi1new);
double pz1new = pt1new/tan(theta1new);
double vz1new = t1->referencePoint().Z();
double vy1new = (vz1new-dz1new)*py1new/pz1new - d01new*px1new/pt1new;
double vx1new = d01new*pt1new/py1new + vy1new*px1new/py1new;
math::XYZPoint refPoint1NEW(vx1new,vy1new,vz1new);
math::XYZVector momentum1NEW(px1new,py1new,pz1new);
const reco::Track *t1_sys = new reco::Track(
t1->chi2(),
t1->ndof(),
refPoint1NEW,
momentum1NEW,
t1->charge(),
t1->covariance(),
t1->algo(),
reco::TrackBase::highPurity,
t1->t0(), t1->beta(), t1->covt0t0(), t1->covBetaBeta()
);
//cout << "ORI 0: " << t0->d0() << " " << t0->dz() << " " << t0->theta() << " " << t0->phi() << " " << t0->pt() << endl;
//cout << "SYS 0: " << t0_sys->d0() << " " << t0_sys->dz() << " " << t0_sys->theta() << " " << t0_sys->phi() << " " << t0_sys->pt() << endl;
//cout << "ORI 1: " << t1->d0() << " " << t1->dz() << " " << t1->theta() << " " << t1->phi() << " " << t1->pt() << endl;
//cout << "SYS 1: " << t1_sys->d0() << " " << t1_sys->dz() << " " << t1_sys->theta() << " " << t1_sys->phi() << " " << t1_sys->pt() << endl;
TLorentzVector muTmp0( lorentzMomentum( myGoodMuonVector[0]->p4() ) );
TLorentzVector muTmp1( lorentzMomentum( myGoodMuonVector[1]->p4() ) );
muon0 = (*theB).build(*t0);
muon1 = (*theB).build(*t1);
muon0_sys = (*theB).build(*t0_sys);
muon1_sys = (*theB).build(*t1_sys);
codeVertexMuo0_ = goodMuons[0];
codeVertexMuo1_ = goodMuons[1];
tks.clear();
tks.push_back(muon0);
tks.push_back(muon1);
tks_sys.clear();
tks_sys.push_back(muon0_sys);
tks_sys.push_back(muon1_sys);
KalmanVertexFitter kalmanMuo(true);
aTransientVertexMuo01 = kalmanMuo.vertex(tks);
KalmanVertexFitter kalmanMuo_sys(true);
aTransientVertexMuo01_sys = kalmanMuo_sys.vertex(tks_sys);
if ( aTransientVertexMuo01.isValid() ) {
//clVertexMuo2_ = TMath::Prob( aTransientVertexMuo01.totalChiSquared(), (int)aTransientVertexMuo01.degreesOfFreedom() );
//cout << aTransientVertexMuo01.totalChiSquared() << " " << aTransientVertexMuo01.degreesOfFreedom() << " " << clVertexMuo2_ << endl;
//cout << ChiSquaredProbability(aTransientVertexMuo01.totalChiSquared(),aTransientVertexMuo01.degreesOfFreedom()) << endl;
// per KF stesso risultato ma per AF dove dof<1 Prob torna sempre 0 e questo ok ma se dof<0 allora risultato NAN
if(aTransientVertexMuo01.degreesOfFreedom()>0) clVertexMuo2_ = ChiSquaredProbability(aTransientVertexMuo01.totalChiSquared(),aTransientVertexMuo01.degreesOfFreedom());
else clVertexMuo2_ = -2;
theSecondaryVertexPositionMuo = aTransientVertexMuo01.position();
//cout << "2 01 " << theSecondaryVertexPositionMuo << endl;
SecVertMuoRec_[0]=aTransientVertexMuo01.position().x();
SecVertMuoRec_[1]=aTransientVertexMuo01.position().y();
SecVertMuoRec_[2]=aTransientVertexMuo01.position().z();
//cout << "Sec Vertex MUO: " << SecVertMuoRec_[0] << " " << SecVertMuoRec_[1] << " " << SecVertMuoRec_[2] << endl;
deltaSecVertPrimVertMuo2_[0] = SecVertMuoRec_[0] - PrimVertRec_[0];
deltaSecVertPrimVertMuo2_[1] = SecVertMuoRec_[1] - PrimVertRec_[1];
deltaSecVertPrimVertMuo2_[2] = SecVertMuoRec_[2] - PrimVertRec_[2];
L_absMuo2_ = sqrt(pow(deltaSecVertPrimVertMuo2_[0],2) + pow(deltaSecVertPrimVertMuo2_[1],2) + pow(deltaSecVertPrimVertMuo2_[2],2));
cov_svMuo2[0][0] = aTransientVertexMuo01.positionError().cxx();
cov_svMuo2[1][0] = aTransientVertexMuo01.positionError().cyx();
cov_svMuo2[2][0] = aTransientVertexMuo01.positionError().czx();
cov_svMuo2[0][1] = cov_svMuo2[1][0];
cov_svMuo2[1][1] = aTransientVertexMuo01.positionError().cyy();
cov_svMuo2[2][1] = aTransientVertexMuo01.positionError().czy();
cov_svMuo2[0][2] = cov_svMuo2[2][0];
cov_svMuo2[1][2] = cov_svMuo2[2][1];
cov_svMuo2[2][2] = aTransientVertexMuo01.positionError().czz();
TLorentzVector muTmp0( lorentzMomentum( myGoodMuonVector[0]->p4() ) );
TLorentzVector muTmp1( lorentzMomentum( myGoodMuonVector[1]->p4() ) );
muon0refitted = aTransientVertexMuo01.refittedTrack(muon0);
muon1refitted = aTransientVertexMuo01.refittedTrack(muon1);
TLorentzVector muRef0( lorentzMomentum( muon0refitted.track().px(), muon0refitted.track().py(), muon0refitted.track().pz(), muMass ) );
TLorentzVector muRef1( lorentzMomentum( muon1refitted.track().px(), muon1refitted.track().py(), muon1refitted.track().pz(), muMass ) );
pxDiMuo = (muRef0+muRef1).Px();
pyDiMuo = (muRef0+muRef1).Py();
pzDiMuo = (muRef0+muRef1).Pz();
piDiMuo = (muRef0+muRef1).P();
//cout << "XTC1: " << pxDiMuo << " " << pyDiMuo << " " << pzDiMuo << " " << piDiMuo << endl;
//cout << "XTC2: " << muTmp0.Pt() << " " << myGoodMuonVector[0]->pt() << endl;
//cout << "XTC3: " << muTmp1.Pt() << " " << myGoodMuonVector[1]->pt() << endl;
//cout << "XTC4: " << myGoodMuonVector.size() << endl;
dimuoMass_ = (muRef0+muRef1).M();
cosphiMuo2_ = (pxDiMuo*deltaSecVertPrimVertMuo2_[0] + pyDiMuo*deltaSecVertPrimVertMuo2_[1] + pzDiMuo*deltaSecVertPrimVertMuo2_[2])/(L_absMuo2_*piDiMuo);
pxDiMuoOri = (muTmp0+muTmp1).Px();
pyDiMuoOri = (muTmp0+muTmp1).Py();
pzDiMuoOri = (muTmp0+muTmp1).Pz();
piDiMuoOri = (muTmp0+muTmp1).P();
//dimuoMassOri_ = (muTmp0+muTmp1).M();
cosphiMuoOri2_ = (pxDiMuoOri*deltaSecVertPrimVertMuo2_[0] + pyDiMuoOri*deltaSecVertPrimVertMuo2_[1] + pzDiMuoOri*deltaSecVertPrimVertMuo2_[2])/(L_absMuo2_*piDiMuoOri);
derivMuo[0] = deltaSecVertPrimVertMuo2_[0]/L_absMuo2_;
derivMuo[1] = deltaSecVertPrimVertMuo2_[1]/L_absMuo2_;
derivMuo[2] = deltaSecVertPrimVertMuo2_[2]/L_absMuo2_;
double LMuo = (pxDiMuo*deltaSecVertPrimVertMuo2_[0] + pyDiMuo*deltaSecVertPrimVertMuo2_[1] + pzDiMuo*deltaSecVertPrimVertMuo2_[2])/piDiMuo;
double sigma_LMuo = 0;
for (int m = 0; m < 3; ++m)
for (int n = 0; n < 3; ++n )
sigma_LMuo += derivMuo[m]*derivMuo[n]*(cov_pv[m][n] + cov_svMuo2[m][n]);
sigma_LMuo = sqrt(sigma_LMuo);
VertexDistance3D vertTool3D;
double distance3D = vertTool3D.distance(primaryVertex, aTransientVertexMuo01).value();
double dist3D_err = vertTool3D.distance(primaryVertex, aTransientVertexMuo01).error();
//cout << "muoVTX prop: " << aTransientVertexMuo01.totalChiSquared() << " " << aTransientVertexMuo01.degreesOfFreedom() << " " << clVertexMuo2_ << " | "
// << distance3D << " " << L_absMuo2_ << " " << LMuo << " | " << dist3D_err << " " << sigma_LMuo << " | " << cosphiMuoOri2_ << endl;
L_3dDistMuo2_ = distance3D;
L_sigmaMuo2_ = LMuo/sigma_LMuo;
}
//cout << "dimuoMassOri : " << dimuoMassOri_ << " dimuoMassRef : " << dimuoMass_ << " cosPhi ori-ref " << cosphiMuoOri2_ << " " << cosphiMuo2_ << endl;
if ( aTransientVertexMuo01_sys.isValid() ) {
//clVertexMuo2_sys_ = TMath::Prob( aTransientVertexMuo01_sys.totalChiSquared(), (int)aTransientVertexMuo01_sys.degreesOfFreedom() );
// per KF stesso risultato ma per AF dove dof<1 Prob torna sempre 0 e questo ok ma se dof<0 allora risultato NAN
if(aTransientVertexMuo01_sys.degreesOfFreedom()>0) clVertexMuo2_sys_ = ChiSquaredProbability(aTransientVertexMuo01_sys.totalChiSquared(),aTransientVertexMuo01_sys.degreesOfFreedom());
else clVertexMuo2_ = -2;
theSecondaryVertexPositionMuo_sys = aTransientVertexMuo01_sys.position();
//cout << "2 01 " << theSecondaryVertexPositionMuo << endl;
SecVertMuoRec_sys_[0]=aTransientVertexMuo01_sys.position().x();
SecVertMuoRec_sys_[1]=aTransientVertexMuo01_sys.position().y();
SecVertMuoRec_sys_[2]=aTransientVertexMuo01_sys.position().z();
//cout << "Sec Vertex MUO: " << SecVertMuoRec_[0] << " " << SecVertMuoRec_[1] << " " << SecVertMuoRec_[2] << endl;
deltaSecVertPrimVertMuo2_sys_[0] = SecVertMuoRec_sys_[0] - PrimVertRec_[0];
deltaSecVertPrimVertMuo2_sys_[1] = SecVertMuoRec_sys_[1] - PrimVertRec_[1];
deltaSecVertPrimVertMuo2_sys_[2] = SecVertMuoRec_sys_[2] - PrimVertRec_[2];
L_absMuo2_sys_ = sqrt(pow(deltaSecVertPrimVertMuo2_sys_[0],2) + pow(deltaSecVertPrimVertMuo2_sys_[1],2) + pow(deltaSecVertPrimVertMuo2_sys_[2],2));
cov_svMuo2_sys[0][0] = aTransientVertexMuo01_sys.positionError().cxx();
cov_svMuo2_sys[1][0] = aTransientVertexMuo01_sys.positionError().cyx();
cov_svMuo2_sys[2][0] = aTransientVertexMuo01_sys.positionError().czx();
cov_svMuo2_sys[0][1] = cov_svMuo2_sys[1][0];
cov_svMuo2_sys[1][1] = aTransientVertexMuo01_sys.positionError().cyy();
cov_svMuo2_sys[2][1] = aTransientVertexMuo01_sys.positionError().czy();
cov_svMuo2_sys[0][2] = cov_svMuo2_sys[2][0];
cov_svMuo2_sys[1][2] = cov_svMuo2_sys[2][1];
cov_svMuo2_sys[2][2] = aTransientVertexMuo01_sys.positionError().czz();
muon0refitted_sys = aTransientVertexMuo01_sys.refittedTrack(muon0_sys);
muon1refitted_sys = aTransientVertexMuo01_sys.refittedTrack(muon1_sys);
TLorentzVector muRef0_sys( lorentzMomentum( muon0refitted_sys.track().px(), muon0refitted_sys.track().py(), muon0refitted_sys.track().pz(), muMass ) );
TLorentzVector muRef1_sys( lorentzMomentum( muon1refitted_sys.track().px(), muon1refitted_sys.track().py(), muon1refitted_sys.track().pz(), muMass ) );
pxDiMuo_sys = (muRef0_sys + muRef1_sys).Px();
pyDiMuo_sys = (muRef0_sys + muRef1_sys).Py();
pzDiMuo_sys = (muRef0_sys + muRef1_sys).Pz();
piDiMuo_sys = (muRef0_sys + muRef1_sys).P();
dimuoMass_sys_ = (muRef0_sys+muRef1_sys).M();
cosphiMuo2_sys_ = (pxDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[0] + pyDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[1] + pzDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[2])/(L_absMuo2_sys_*piDiMuo_sys);
derivMuo_sys[0] = deltaSecVertPrimVertMuo2_sys_[0]/L_absMuo2_sys_;
derivMuo_sys[1] = deltaSecVertPrimVertMuo2_sys_[1]/L_absMuo2_sys_;
derivMuo_sys[2] = deltaSecVertPrimVertMuo2_sys_[2]/L_absMuo2_sys_;
double LMuo_sys = (pxDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[0] + pyDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[1] + pzDiMuo_sys*deltaSecVertPrimVertMuo2_sys_[2])/piDiMuo_sys;
double sigma_LMuo_sys = 0;
for (int m = 0; m < 3; ++m)
for (int n = 0; n < 3; ++n )
sigma_LMuo_sys += derivMuo_sys[m]*derivMuo_sys[n]*(cov_pv[m][n] + cov_svMuo2_sys[m][n]);
sigma_LMuo_sys = sqrt(sigma_LMuo_sys);
VertexDistance3D vertTool3D_sys;
double distance3D_sys = vertTool3D_sys.distance(primaryVertex, aTransientVertexMuo01_sys).value();
//double dist3D_err = vertTool3D.distance(primaryVertex, aTransientVertexMuo01).error();
//cout << "muoVTX prop: " << aTransientVertexMuo01.totalChiSquared() << " " << aTransientVertexMuo01.degreesOfFreedom() << " " << clVertexMuo2_ << " | "
// << distance3D << " " << L_absMuo2_ << " " << LMuo << " | " << dist3D_err << " " << sigma_LMuo << " | " << cosphiMuoOri2_ << endl;
L_3dDistMuo2_sys_ = distance3D_sys;
L_sigmaMuo2_sys_ = LMuo_sys/sigma_LMuo_sys;
}
//cout << "NO SYS: " << SecVertMuoRec_[0] << ":" << SecVertMuoRec_[1] << ":" << SecVertMuoRec_[2] << " " << cosphiMuo2_ << " " << L_sigmaMuo2_ << endl;
//cout << "SI SYS: " << SecVertMuoRec_sys_[0] << ":" << SecVertMuoRec_sys_[1] << ":" << SecVertMuoRec_sys_[2] << " " << cosphiMuo2_sys_ << " " << L_sigmaMuo2_sys_ << endl;
delete t0_sys;
delete t1_sys;
}
if( clVertexMuo2_ > 0 ) {
//cout << "------- MUO SV: " << clVertexMuo2_ << " " << theSecondaryVertexPositionMuo << endl;
TrajectoryStateClosestToPoint pcaState0 = muon0refitted.trajectoryStateClosestToPoint(theSecondaryVertexPositionMuo);
const GlobalPoint& pcaPosition0 = pcaState0.position();
//cout << "TCP0 " << pcaPosition0 << endl;
auto deltaSecVertexPCA0 = pcaPosition0 - theSecondaryVertexPositionMuo;
muo0dis_ = sqrt(deltaSecVertexPCA0*deltaSecVertexPCA0);
//cout << "D0 " << deltaSecVertexPCA0 << " " << muo0dis_ << endl;
TrajectoryStateClosestToPoint pcaState1 = muon1refitted.trajectoryStateClosestToPoint(theSecondaryVertexPositionMuo);
const GlobalPoint& pcaPosition1 = pcaState1.position();
//cout << "TCP1 " << pcaPosition1 << endl;
auto deltaSecVertexPCA1 = pcaPosition1 - theSecondaryVertexPositionMuo;
muo1dis_ = sqrt(deltaSecVertexPCA1*deltaSecVertexPCA1);
//cout << "D1 " << deltaSecVertexPCA1 << " " << muo1dis_ << endl;
}
//-------------------------//
fTree->Fill();
gen->Clear("C");
muo->Clear("C");
muotrack->Clear("C");
acceptedTriggers->Clear("C");
}
MakeFinalMuoTreeMINIAOD::MakeFinalMuoTreeMINIAOD(const edm::ParameterSet& cfg):
genEvtInfoToken_ (consumes<GenEventInfoProduct> (cfg.getParameter<edm::InputTag>("genEvtInfo"))),
genParticleToken_ (consumes<std::vector<pat::PackedGenParticle> > (cfg.getParameter<edm::InputTag>("GenParticles"))),
mcTruthCollectionToken_ (consumes<edm::View<reco::GenParticle> > (cfg.getParameter<edm::InputTag>("mcTruthCollection"))),
puSummaryToken_ (consumes<std::vector<PileupSummaryInfo> > (cfg.getParameter<edm::InputTag>("pileupSummaryInfos"))),
vtxToken_ (consumes<reco::VertexCollection> (cfg.getParameter<edm::InputTag>("vertices"))),
metToken_ (consumes<std::vector<pat::MET> > (cfg.getParameter<edm::InputTag>("mets"))),
muonToken_ (consumes<std::vector<pat::Muon> > (cfg.getParameter<edm::InputTag>("muons"))),
metFiltersBitToken_ (consumes<edm::TriggerResults> (cfg.getParameter<edm::InputTag>("metfilters"))),
triggerBitToken_ (consumes<edm::TriggerResults> (cfg.getParameter<edm::InputTag>("bits"))),
pfCandToken_ (consumes<std::vector<pat::PackedCandidate> > (cfg.getParameter<edm::InputTag>("pfcands"))),
dedxMapToken_ (consumes<reco::DeDxHitInfoAss> (cfg.getParameter<edm::InputTag>("DeDxHitMap"))),
dedxHitInfoToken_ (consumes<std::vector<reco::DeDxHitInfo> > (cfg.getParameter<edm::InputTag>("DeDxHitInfo"))),
beamSpotToken_ (consumes<reco::BeamSpot> (cfg.getParameter<edm::InputTag>("offlineBeamSpot"))),
transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
hltPrescaleProvider_ (cfg, consumesCollector(), *this)
{
prefweight_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProb"));
prefweightup_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbUp"));
prefweightdown_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbDown"));
prefweightECAL_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbECAL"));
prefweightECALup_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbECALUp"));
prefweightECALdown_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbECALDown"));
prefweightMUON_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbMuon"));
prefweightMUONup_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbMuonUp"));
prefweightMUONdown_token = consumes< double >(edm::InputTag("prefiringweight:nonPrefiringProbMuonDown"));
BadPFMuonDzFilterToken_ = consumes< bool >(edm::InputTag("BadPFMuonFilterUpdateDz"));
sysRun_ = cfg.getParameter<int>("sysRun");
era = cfg.getParameter<std::string>("era");
theOutFileName = cfg.getParameter<std::string>("OutputFileName");
outFile = new TFile(theOutFileName.c_str(), "RECREATE", "");
outFile->cd();
fTree = new TTree ("ana", "CMSSW anaFind tree");
fTree->SetAutoSave(500000); // Save every ~1000 events
//fTree->SetMaxTreeSize(20000000000);
fTree->Branch("event", &eventNb, "event/i");
fTree->Branch("run", &runNb, "run/i");
fTree->Branch("lumi", &lumiBlock, "lumi/i");
fTree->Branch("nvtx", &nvtx, "nvtx/I");
fTree->Branch("nputrue", &nputrue, "nputrue/F");
fTree->Branch("genweight", &genweight, "genweight/F");
fTree->Branch("metFilters", &metFilters, "metFilters/O");
fTree->Branch("prefiringweightECAL", &prefiringweightECAL_, "prefiringweightECAL/D");
fTree->Branch("prefiringweightECALup", &prefiringweightECALup_, "prefiringweightECALup/D");
fTree->Branch("prefiringweightECALdown", &prefiringweightECALdown_, "prefiringweightECALdown/D");
fTree->Branch("prefiringweightMUON", &prefiringweightMUON_, "prefiringweightMUON/D");
fTree->Branch("prefiringweightMUONup", &prefiringweightMUONup_, "prefiringweightMUONup/D");
fTree->Branch("prefiringweightMUONdown", &prefiringweightMUONdown_, "prefiringweightMUONdown/D");
fTree->Branch("metx", &pfMETx, "metx/F");
fTree->Branch("mety", &pfMETy, "mety/F");
fTree->Branch("met", &pfMET, "met/F");
fTree->Branch("ht", &ht, "ht/F");
acceptedTriggers = new TClonesArray("TObjString", 10000);
gen = new TClonesArray("TLorentzVector", 10000);
muo = new TClonesArray("TLorentzVector", 10000);
muotrack = new TClonesArray("TLorentzVector", 10000);
fTree->Branch("nGen", &nGen_, "nGen/I");
fTree->Branch("idGen", &idGen_, "idGen[nGen]/I");
fTree->Branch("PrimVertGen", &PrimVertGen_, "PrimVertGen[3]/D");
fTree->Branch("PrimVertRec", &PrimVertRec_, "PrimVertRec[3]/D");
fTree->Branch("PrimVertRefRec", &PrimVertRefRec_, "PrimVertRefRec[3]/D");
fTree->Branch("SecVertGen", &SecVertGen_, "SecVertGen[3]/D");
fTree->Branch("SecVertMuoRec", &SecVertMuoRec_, "SecVertMuoRec[3]/D");
fTree->Branch("iAcceptedTriggers", &iAcceptedTriggers_, "iAcceptedTriggers/I");
fTree->Branch("nMuo", &nMuo_, "nMuo/I");
fTree->Branch("idMuo", &idMuo_, "idMuo[nMuo]/I");
fTree->Branch("nlMuo", &nlMuo_, "nlMuo[nMuo]/I");
fTree->Branch("fHLTMuo", &fHLTMuo_, "fHLTMuo[nMuo]/I");
fTree->Branch("sip3dMuo", &sip3dMuo_, "sip3dMuo[nMuo]/F");
fTree->Branch("isoMuoPF", &isoMuoPF_, "isoMuoPF[nMuo]/F");
fTree->Branch("isoMuoTR", &isoMuoTR_, "isoMuoTR[nMuo]/F");
fTree->Branch("dxyMuo", &dxyMuo_, "dxyMuo[nMuo]/F");
fTree->Branch("dzMuo", &dzMuo_, "dzMuo[nMuo]/F");
fTree->Branch("dedxMuo", &dedxMuo_, "dedxMuo[nMuo]/F");
fTree->Branch("dedxMuoNhits", &dedxMuoNhits_, "dedxMuoNhits[nMuo]/I");
fTree->Branch("ptErrorMuo", &ptErrorMuo_, "ptErrorMuo[nMuo]/F");
fTree->Branch("clVertexMuo2", &clVertexMuo2_, "clVertexMuo2/D");
fTree->Branch("codeVertexMuo0", &codeVertexMuo0_, "codeVertexMuo0/I");
fTree->Branch("codeVertexMuo1", &codeVertexMuo1_, "codeVertexMuo1/I");
fTree->Branch("L_sigmaMuo2", &L_sigmaMuo2_, "L_sigmaMuo2/D");
fTree->Branch("L_absMuo2", &L_absMuo2_, "L_absMuo2/D");
fTree->Branch("L_3dDistMuo2", &L_3dDistMuo2_, "L_3dDistMuo2/D");
fTree->Branch("cosphiMuo2", &cosphiMuo2_, "cosphiMuo2/D");
fTree->Branch("cosphiMuoOri2", &cosphiMuoOri2_, "cosphiMuoOri2/D");
fTree->Branch("dimuoMass", &dimuoMass_, "dimuoMass/F");
fTree->Branch("dimuoCharge", &dimuoCharge_, "dimuoCharge/I");
fTree->Branch("deltaSecVertPrimVertMuo2", &deltaSecVertPrimVertMuo2_, "deltaSecVertPrimVertMuo2[3]/D");
fTree->Branch("muo0dis", &muo0dis_, "muo0dis/D");
fTree->Branch("muo1dis", &muo1dis_, "muo1dis/D");
fTree->Branch("clVertexMuo2_sys", &clVertexMuo2_sys_, "clVertexMuo2_sys/D");
fTree->Branch("L_sigmaMuo2_sys", &L_sigmaMuo2_sys_, "L_sigmaMuo2_sys/D");
fTree->Branch("L_absMuo2_sys", &L_absMuo2_sys_, "L_absMuo2_sys/D");
fTree->Branch("L_3dDistMuo2_sys", &L_3dDistMuo2_sys_, "L_3dDistMuo2_sys/D");
fTree->Branch("cosphiMuo2_sys", &cosphiMuo2_sys_, "cosphiMuo2_sys/D");
fTree->Branch("dimuoMass_sys", &dimuoMass_sys_, "dimuoMass_sys/F");
fTree->Branch("deltaSecVertPrimVertMuo2_sys", &deltaSecVertPrimVertMuo2_sys_, "deltaSecVertPrimVertMuo2_sys[3]/D");
fTree->Branch("acceptedTriggers", "TClonesArray", &acceptedTriggers, 32000, 0);
fTree->Branch("gen", "TClonesArray", &gen, 32000, 0);
fTree->Branch("muo", "TClonesArray", &muo, 32000, 0);
fTree->Branch("muotrack", "TClonesArray", &muotrack, 32000, 0);
}
///////////////////////////////////////////////////////////////
// called at end
///////////////////////////////////////////////////////////////
void MakeFinalMuoTreeMINIAOD::endJob() {
TFile *actualFile = fTree->GetCurrentFile();
actualFile->Write();
actualFile->Close();
}
///////////////////////////////////////////////////////////////////////////////
void MakeFinalMuoTreeMINIAOD::beginRun(const edm::Run &iRun, const edm::EventSetup &iSetup) {
//std::cout << "calling init(" << "iRun" << ", " << "iSetup" << ", " << ", "<< "changed_" << ") in beginRun()" << std::endl;
bool changed(true);
if (!hltPrescaleProvider_.init(iRun, iSetup, "HLT", changed)) {
edm::LogError("CandidateTriggerObjectProducer") << "Error! Can't initialize HLTConfigProvider";
throw cms::Exception("HLTConfigProvider::init() returned non 0");
}
}
///////////////////////////////////////////////////////////////////////////////
// Returns Lorentz-vector of PF particle
//////////////////////////////////////////////////////////////////////////////
TLorentzVector MakeFinalMuoTreeMINIAOD::lorentzMomentum(const LorentzVector pfcl) const {
double ereco = pfcl.energy();
double pxreco = pfcl.px();
double pyreco = pfcl.py();
double pzreco = pfcl.pz();
// cout << "pfcl = " << pxreco << " " << pyreco << " " << pzreco << " " << ereco << endl;
TLorentzVector lrzpreco(pxreco, pyreco, pzreco, ereco);
return lrzpreco;
}
///////////////////////////////////////////////////////////////////////////////
// Returns Lorentz-vector of PF particle
//////////////////////////////////////////////////////////////////////////////
TLorentzVector MakeFinalMuoTreeMINIAOD::lorentzMomentum(double pxreco, double pyreco, double pzreco, double mass) const {
double ereco = sqrt(pxreco*pxreco+pyreco*pyreco+pzreco*pzreco+mass*mass);
TLorentzVector lrzpreco(pxreco, pyreco, pzreco, ereco);
return lrzpreco;
}
// copiato da https://github.com/cms-sw/cmssw/blob/23be6cb426861c3cc2b1b9f58c32ee679839bc0a/DQM/TrackingMonitor/src/VertexMonitor.cc#L376
void MakeFinalMuoTreeMINIAOD::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
//The following says we do not know what parameters are allowed so do no validation
// Please change this to state exactly what you do use, even if it is no parameters
edm::ParameterSetDescription desc;
desc.add<unsigned int>("stageL1Trigger", 2);
desc.add<edm::InputTag>("genEvtInfo",edm::InputTag("generator"));
desc.add<edm::InputTag>("DeDxHitInfo",edm::InputTag("isolatedTracks"));
desc.add<edm::InputTag>("DeDxHitMap",edm::InputTag("isolatedTracks"));
desc.add<edm::InputTag>("GenParticles",edm::InputTag("packedGenParticles"));
desc.add<string>("OutputFileName","myTree.root");
desc.add<edm::InputTag>("bits",edm::InputTag("TriggerResults","","HLT"));
desc.add<string>("era","2022b");
desc.add<edm::InputTag>("mcTruthCollection",edm::InputTag("prunedGenParticles"));
desc.add<edm::InputTag>("metfilters",edm::InputTag("TriggerResults","","RECO"));
desc.add<edm::InputTag>("mets",edm::InputTag("slimmedMETs"));
desc.add<edm::InputTag>("muons",edm::InputTag("slimmedMuons"));
desc.add<edm::InputTag>("offlineBeamSpot",edm::InputTag("offlineBeamSpot"));
desc.add<edm::InputTag>("pfcands",edm::InputTag("packedPFCandidates"));
desc.add<edm::InputTag>("pileupSummaryInfos",edm::InputTag("slimmedAddPileupInfo"));
desc.add<int>("sysRun",0);
desc.add<edm::InputTag>("vertices",edm::InputTag("offlineSlimmedPrimaryVertices"));
descriptions.add("MakeFinalMuoTreeMINIAOD", desc);
}
//--------------------------------------------------------------------------------------------------------------
//define this as a plug-in
DEFINE_FWK_MODULE(MakeFinalMuoTreeMINIAOD);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment