MomentumScaleCalibration/plugins/MuScleFitFilter.cc

0001 // #define DEBUG
0002 // System include files
0003 // --------------------
0004
0005 // User include files
0006 // ------------------
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDFilter.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/MuonReco/interface/Muon.h"
0015 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0016 #include "DataFormats/Candidate/interface/Candidate.h"
0017 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0018 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0019
0020 #include <CLHEP/Vector/LorentzVector.h>
0021
0022 // For file output
0023 // ---------------
0024 #include <fstream>
0025 #include <sstream>
0026 #include <cmath>
0027 #include <memory>
0028
0029 #include <vector>
0030
0031 #include "TRandom.h"
0032
0033 // Class declaration
0034 // -----------------
0035
0036 class MuScleFitFilter : public edm::stream::EDFilter<> {
0037 public:
0038   explicit MuScleFitFilter(const edm::ParameterSet&);
0039   ~MuScleFitFilter() override;
0040
0041 private:
0042   bool filter(edm::Event&, const edm::EventSetup&) override;
0043
0044   // Member data
0045   // -----------
0046   int eventsRead;
0047   int eventsWritten;
0048   bool debug;
0049   int theMuonType;
0050   std::vector<double> Mmin;
0051   std::vector<double> Mmax;
0052   int maxWrite;
0053   unsigned int minimumMuonsNumber;
0054
0055   // Collections labels
0056   // ------------------
0057   edm::InputTag theMuonLabel;
0058   edm::EDGetTokenT<reco::MuonCollection> theGlbMuonsToken;
0059   edm::EDGetTokenT<reco::TrackCollection> theSaMuonsToken;
0060   edm::EDGetTokenT<reco::TrackCollection> theTracksToken;
0061 };
0062
0063 // Static data member definitions
0064 // ------------------------------
0065 const double Mmu2 = 0.011163612;  // Squared muon mass
0066
0067 // Constructor
0068 // -----------
0069 MuScleFitFilter::MuScleFitFilter(const edm::ParameterSet& iConfig) {
0070   debug = iConfig.getUntrackedParameter<bool>("debug", false);
0071
0072   if (debug)
0073     edm::LogPrint("MuScleFitFilter") << "Constructor" << std::endl;
0074
0075   // Parameters
0076   // ----------
0077   //ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0078   //theService = new MuonServiceProxy(serviceParameters);
0079   theMuonLabel = iConfig.getParameter<edm::InputTag>("MuonLabel");
0080   theGlbMuonsToken = mayConsume<reco::MuonCollection>(theMuonLabel);
0081   theSaMuonsToken = mayConsume<reco::TrackCollection>(theMuonLabel);
0082   theTracksToken = mayConsume<reco::TrackCollection>(theMuonLabel);
0083   theMuonType = iConfig.getParameter<int>("muonType");
0084
0085   Mmin = iConfig.getUntrackedParameter<std::vector<double> >("Mmin");
0086   Mmax = iConfig.getUntrackedParameter<std::vector<double> >("Mmax");
0087   maxWrite = iConfig.getUntrackedParameter<int>("maxWrite", 100000);
0088
0089   minimumMuonsNumber = iConfig.getUntrackedParameter<unsigned int>("minimumMuonsNumber", 2);
0090
0091   // The must have the same size and they must not be empty, otherwise abort.
0092   if (!(Mmin.size() == Mmax.size() && !Mmin.empty()))
0093     abort();
0094
0095   // Count the total number of analyzed and written events
0096   // -----------------------------------------------------
0097   eventsRead = 0;
0098   eventsWritten = 0;
0099 }
0100
0101 // Destructor
0102 // ----------
0103 MuScleFitFilter::~MuScleFitFilter() {
0104   edm::LogPrint("MuScleFitFilter") << "Total number of events read    = " << eventsRead << std::endl;
0105   edm::LogPrint("MuScleFitFilter") << "Total number of events written = " << eventsWritten << std::endl;
0106 }
0107
0108 // Member functions
0109 // ----------------
0110
0111 // Method called for each event
0112 // ----------------------------
0113 bool MuScleFitFilter::filter(edm::Event& event, const edm::EventSetup& iSetup) {
0114   // Cut the crap if we have stored enough stuff
0115   // -------------------------------------------
0116   if (maxWrite != -1 && eventsWritten >= maxWrite)
0117     return false;
0118
0119   // Get the RecTrack and the RecMuon collection from the event
0120   // ----------------------------------------------------------
0121   std::unique_ptr<reco::MuonCollection> muons(new reco::MuonCollection());
0122
0123   if (debug)
0124     edm::LogPrint("MuScleFitFilter") << "Looking for muons of the right kind" << std::endl;
0125
0126   if (theMuonType == 1) {  // GlobalMuons
0127
0128     // Global muons:
0129     // -------------
0130     edm::Handle<reco::MuonCollection> glbMuons;
0131     if (debug)
0132       edm::LogPrint("MuScleFitFilter") << "Handle defined" << std::endl;
0133     event.getByToken(theGlbMuonsToken, glbMuons);
0134     if (debug)
0135       edm::LogPrint("MuScleFitFilter") << "Global muons: " << glbMuons->size() << std::endl;
0136
0137     // Store the muon
0138     // --------------
0139     reco::MuonCollection::const_iterator glbMuon;
0140     for (glbMuon = glbMuons->begin(); glbMuon != glbMuons->end(); ++glbMuon) {
0141       muons->push_back(*glbMuon);
0142       if (debug) {
0143         edm::LogPrint("MuScleFitFilter") << "  Reconstructed muon: pT = " << glbMuon->p4().Pt()
0144                                          << "  Eta = " << glbMuon->p4().Eta() << std::endl;
0145       }
0146     }
0147   } else if (theMuonType == 2) {  // StandaloneMuons
0148
0149     // Standalone muons:
0150     // -----------------
0151     edm::Handle<reco::TrackCollection> saMuons;
0152     event.getByToken(theSaMuonsToken, saMuons);
0153     if (debug)
0154       edm::LogPrint("MuScleFitFilter") << "Standalone muons: " << saMuons->size() << std::endl;
0155
0156     // Store the muon
0157     // --------------
0158     reco::TrackCollection::const_iterator saMuon;
0159     for (saMuon = saMuons->begin(); saMuon != saMuons->end(); ++saMuon) {
0160       reco::Muon muon;
0161       double energy = sqrt(saMuon->p() * saMuon->p() + Mmu2);
0162       math::XYZTLorentzVector p4(saMuon->px(), saMuon->py(), saMuon->pz(), energy);
0163       muon.setP4(p4);
0164       muon.setCharge(saMuon->charge());
0165       muons->push_back(muon);
0166     }
0167   } else if (theMuonType == 3) {  // Tracker tracks
0168
0169     // Tracks:
0170     // -------
0171     edm::Handle<reco::TrackCollection> tracks;
0172     event.getByToken(theTracksToken, tracks);
0173     if (debug)
0174       edm::LogPrint("MuScleFitFilter") << "Tracker tracks: " << tracks->size() << std::endl;
0175
0176     // Store the muon
0177     // -------------
0178     reco::TrackCollection::const_iterator track;
0179     for (track = tracks->begin(); track != tracks->end(); ++track) {
0180       reco::Muon muon;
0181       double energy = sqrt(track->p() * track->p() + Mmu2);
0182       math::XYZTLorentzVector p4(track->px(), track->py(), track->pz(), energy);
0183       muon.setP4(p4);
0184       muon.setCharge(track->charge());
0185       muons->push_back(muon);
0186     }
0187   } else {
0188     edm::LogPrint("MuScleFitFilter") << "Wrong muon type! Aborting." << std::endl;
0189     abort();
0190   }
0191
0192   // Loop on RecMuon and reconstruct the resonance
0193   // ---------------------------------------------
0194   reco::MuonCollection::const_iterator muon1;
0195   reco::MuonCollection::const_iterator muon2;
0196
0197   bool resfound = false;
0198
0199   // Require at least N muons of the selected type.
0200   if (muons->size() >= minimumMuonsNumber) {
0201     for (muon1 = muons->begin(); muon1 != muons->end(); ++muon1) {
0202       if (debug) {
0203         edm::LogPrint("MuScleFitFilter") << "  Reconstructed muon: pT = " << muon1->p4().Pt()
0204                                          << "  Eta = " << muon1->p4().Eta() << std::endl;
0205       }
0206
0207       // Recombine all the possible Z from reconstructed muons
0208       // -----------------------------------------------------
0209       if (muons->size() > 1) {
0210         for (muon2 = muon1 + 1; muon2 != muons->end(); ++muon2) {
0211           if (((*muon1).charge() * (*muon2).charge()) < 0) {  // This also gets rid of muon1==muon2
0212             //    reco::Particle::LorentzVector Z (muonCorr1 + muonCorr2);
0213             reco::Particle::LorentzVector Z(muon1->p4() + muon2->p4());
0214             // Loop on all the cuts on invariant mass.
0215             // If it passes at least one of the cuts, the event will be accepted.
0216             // ------------------------------------------------------------------
0217             std::vector<double>::const_iterator mMinCut = Mmin.begin();
0218             std::vector<double>::const_iterator mMaxCut = Mmax.begin();
0219             for (; mMinCut != Mmin.end(); ++mMinCut, ++mMaxCut) {
0220               // When the two borders are -1 do not cut.
0221               if (*mMinCut == *mMaxCut && *mMaxCut == -1) {
0222                 resfound = true;
0223                 if (debug) {
0224                   edm::LogPrint("MuScleFitFilter")
0225                       << "Acceptiong event because mMinCut = " << *mMinCut << " = mMaxCut = " << *mMaxCut << std::endl;
0226                 }
0227               } else if (Z.mass() > *mMinCut && Z.mass() < *mMaxCut) {
0228                 resfound = true;
0229                 if (debug) {
0230                   edm::LogPrint("MuScleFitFilter") << "One particle found with mass = " << Z.mass() << std::endl;
0231                 }
0232               }
0233             }
0234           }
0235         }
0236       } else if (debug) {
0237         edm::LogPrint("MuScleFitFilter") << "Not enough reconstructed muons to make a resonance" << std::endl;
0238       }
0239     }
0240   } else if (debug) {
0241     edm::LogPrint("MuScleFitFilter") << "Skipping event because muons = " << muons->size() << " < "
0242                                      << "minimumMuonsNumber(" << minimumMuonsNumber << ")" << std::endl;
0243   }
0244
0245   // Store the event if it has a dimuon pair with mass within defined boundaries
0246   // ---------------------------------------------------------------------------
0247   bool write = false;
0248   eventsRead++;
0249   if (resfound) {
0250     write = true;
0251     eventsWritten++;
0252   }
0253   return write;
0254 }
0255
0256 #include "FWCore/Framework/interface/MakerMacros.h"
0257 DEFINE_FWK_MODULE(MuScleFitFilter);