Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-07 22:41:46

0001 
0002 // -*- C++ -*-
0003 //
0004 // Package:    Configuration/Skimming
0005 // Class:      LeptonSkimming
0006 //
0007 /**\class LeptonSkimming LeptonSkimming.cc Configuration/Skimming/src/LeptonSkimming.cc
0008 
0009  Description: [one line class summary]
0010 
0011  Implementation:
0012      [Notes on implementation]
0013 */
0014 //
0015 // Original Author:  Georgios Karathanasis georgios.karathanasis@cern.ch
0016 //         Created:  Thu, 29 Nov 2018 15:23:09 GMT
0017 //
0018 //
0019 
0020 #include "Configuration/Skimming/interface/LeptonSkimming.h"
0021 
0022 using namespace edm;
0023 using namespace reco;
0024 using namespace std;
0025 
0026 LeptonSkimming::LeptonSkimming(const edm::ParameterSet& iConfig)
0027     : electronsToken_(consumes<std::vector<reco::GsfElectron>>(iConfig.getParameter<edm::InputTag>("electrons"))),
0028       eleBWPToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("eleBiasedWP"))),
0029       eleUnBWPToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("eleUnbiasedWP"))),
0030       muonsToken_(consumes<std::vector<reco::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))),
0031       Tracks_(consumes<std::vector<reco::Track>>(iConfig.getParameter<edm::InputTag>("tracks"))),
0032       vtxToken_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("vertices"))),
0033       beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0034       conversionsToken_(consumes<reco::ConversionCollection>(iConfig.getParameter<edm::InputTag>("conversions"))),
0035 
0036       trgresultsToken_(consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("triggerresults"))),
0037       trigobjectsToken_(consumes<trigger::TriggerEvent>(iConfig.getParameter<edm::InputTag>("triggerobjects"))),
0038       bFieldToken_(esConsumes()),
0039       HLTFilter_(iConfig.getParameter<vector<string>>("HLTFilter")),
0040       HLTPath_(iConfig.getParameter<vector<string>>("HLTPath")) {
0041   edm::ParameterSet runParameters = iConfig.getParameter<edm::ParameterSet>("RunParameters");
0042   //mu matching with trigger
0043   MuTrgMatchCone = runParameters.getParameter<double>("MuTrgMatchCone");
0044   SkipIfNoMuMatch = runParameters.getParameter<bool>("SkipIfNoMuMatch");
0045   //track cuts
0046   PtTrack_Cut = runParameters.getParameter<double>("PtTrack_Cut");
0047   EtaTrack_Cut = runParameters.getParameter<double>("EtaTrack_Cut");
0048   MinChi2Track_Cut = runParameters.getParameter<double>("MinChi2Track_Cut");
0049   MaxChi2Track_Cut = runParameters.getParameter<double>("MaxChi2Track_Cut");
0050   MuTrkMinDR_Cut = runParameters.getParameter<double>("MuTrkMinDR_Cut");
0051   TrackMuDz_Cut = runParameters.getParameter<double>("TrackMuDz_Cut");
0052   TrgExclusionCone = runParameters.getParameter<double>("TrgExclusionCone");
0053   //lepton cuts
0054   PtMu_Cut = runParameters.getParameter<double>("PtMu_Cut");
0055   PtEl_Cut = runParameters.getParameter<double>("PtEl_Cut");
0056   QualMu_Cut = runParameters.getParameter<double>("QualMu_Cut");
0057   MuTrgExclusionCone = runParameters.getParameter<double>("MuTrgExclusionCone");
0058   ElTrgExclusionCone = runParameters.getParameter<double>("ElTrgExclusionCone");
0059   TrkObjExclusionCone = runParameters.getParameter<double>("TrkObjExclusionCone");
0060   MuTrgMuDz_Cut = runParameters.getParameter<double>("MuTrgMuDz_Cut");
0061   ElTrgMuDz_Cut = runParameters.getParameter<double>("ElTrgMuDz_Cut");
0062   BiasedWP = runParameters.getParameter<double>("BiasedWP");
0063   UnbiasedWP = runParameters.getParameter<double>("UnbiasedWP");
0064   SkimOnlyMuons = runParameters.getParameter<bool>("SkimOnlyMuons");
0065   SkimOnlyElectrons = runParameters.getParameter<bool>("SkimOnlyElectrons");
0066   //pair cuts
0067   MaxMee_Cut = runParameters.getParameter<double>("MaxMee_Cut");
0068   MinMee_Cut = runParameters.getParameter<double>("MinMee_Cut");
0069   Probee_Cut = runParameters.getParameter<double>("Probee_Cut");
0070   Cosee_Cut = runParameters.getParameter<double>("Cosee_Cut");
0071   EpairZvtx_Cut = runParameters.getParameter<double>("EpairZvtx_Cut");
0072   //kaon
0073   PtKTrack_Cut = runParameters.getParameter<double>("PtKTrack_Cut");
0074   TrackSdxy_Cut = runParameters.getParameter<double>("TrackSdxy_Cut");
0075   Ksdxy_Cut = runParameters.getParameter<double>("Ksdxy_Cut");
0076   //triplet
0077   MaxMB_Cut = runParameters.getParameter<double>("MaxMB_Cut");
0078   MinMB_Cut = runParameters.getParameter<double>("MinMB_Cut");
0079   ProbeeK_Cut = runParameters.getParameter<double>("ProbeeK_Cut");
0080   CoseeK_Cut = runParameters.getParameter<double>("CoseeK_Cut");
0081   SLxy_Cut = runParameters.getParameter<double>("SLxy_Cut");
0082   PtB_Cut = runParameters.getParameter<double>("PtB_Cut");
0083 
0084   ObjPtLargerThanTrack = runParameters.getParameter<bool>("ObjPtLargerThanTrack");
0085 }
0086 
0087 LeptonSkimming::~LeptonSkimming() {
0088   // do anything here that needs to be done at destruction time
0089   // (e.g. close files, deallocate resources etc.)
0090 }
0091 
0092 bool LeptonSkimming::hltFired(const edm::Event& iEvent, const edm::EventSetup& iSetup, std::vector<string> HLTPath) {
0093   using namespace std;
0094   using namespace edm;
0095   using namespace reco;
0096   using namespace trigger;
0097 
0098   edm::Handle<edm::TriggerResults> trigResults;
0099   iEvent.getByToken(trgresultsToken_, trigResults);
0100   bool fire = false;
0101   if (trigResults.failedToGet())
0102     return false;
0103   for (unsigned int ip = 0; ip < HLTPath.size(); ip++) {
0104     int N_Triggers = trigResults->size();
0105     const edm::TriggerNames& trigName = iEvent.triggerNames(*trigResults);
0106     for (int i_Trig = 0; i_Trig < N_Triggers; ++i_Trig) {
0107       if (!trigResults->accept(i_Trig))
0108         continue;
0109       const std::string& TrigPath = trigName.triggerName(i_Trig);
0110       if (TrigPath.find(HLTPath[ip]) != std::string::npos)
0111         fire = true;
0112     }
0113   }
0114   return fire;
0115 }
0116 
0117 std::array<float, 5> LeptonSkimming::hltObject(const edm::Event& iEvent,
0118                                                const edm::EventSetup& iSetup,
0119                                                std::vector<string> Seed) {
0120   using namespace std;
0121   using namespace edm;
0122   using namespace reco;
0123   using namespace trigger;
0124 
0125   edm::Handle<trigger::TriggerEvent> triggerObjectsSummary;
0126   iEvent.getByToken(trigobjectsToken_, triggerObjectsSummary);
0127   trigger::TriggerObjectCollection selectedObjects;
0128 
0129   std::vector<std::array<float, 5>> max_per_trigger;
0130 
0131   for (unsigned int ipath = 0; ipath < Seed.size(); ipath++) {
0132     std::vector<std::array<float, 5>> tot_tr_obj_pt_eta_phi;
0133     if (!triggerObjectsSummary.isValid())
0134       continue;
0135     size_t filterIndex = (*triggerObjectsSummary).filterIndex(InputTag(Seed[ipath], "", "HLT"));
0136     trigger::TriggerObjectCollection allTriggerObjects = triggerObjectsSummary->getObjects();
0137     if (filterIndex < (*triggerObjectsSummary).sizeFilters()) {
0138       const trigger::Keys& keys = (*triggerObjectsSummary).filterKeys(filterIndex);
0139       for (size_t j = 0; j < keys.size(); j++) {
0140         const trigger::TriggerObject& foundObject = (allTriggerObjects)[keys[j]];
0141         std::array<float, 5> tr_obj_pt_eta_phi;
0142         if (fabs(foundObject.id()) != 13)
0143           continue;
0144         tr_obj_pt_eta_phi[0] = foundObject.pt();
0145         tr_obj_pt_eta_phi[1] = foundObject.eta();
0146         tr_obj_pt_eta_phi[2] = foundObject.phi();
0147         tr_obj_pt_eta_phi[3] = foundObject.id() / fabs(foundObject.id());
0148         tot_tr_obj_pt_eta_phi.push_back(tr_obj_pt_eta_phi);
0149       }
0150     }
0151     //take the max per hlt
0152     if (!tot_tr_obj_pt_eta_phi.empty()) {
0153       std::sort(tot_tr_obj_pt_eta_phi.begin(),
0154                 tot_tr_obj_pt_eta_phi.end(),
0155                 [](const std::array<float, 5>& a, const std::array<float, 5>& b) { return a[0] > b[0]; });
0156       max_per_trigger.push_back(tot_tr_obj_pt_eta_phi.at(0));
0157     }
0158   }
0159   //we know that at least a trigger fired
0160   //find the total max
0161   std::sort(max_per_trigger.begin(),
0162             max_per_trigger.end(),
0163             [](const std::array<float, 5>& a, const std::array<float, 5>& b) { return a[0] > b[0]; });
0164   return max_per_trigger.at(0);
0165 }
0166 
0167 // ------------ method called on each new Event  ------------
0168 bool LeptonSkimming::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0169   using namespace std;
0170   using namespace edm;
0171   using namespace reco;
0172   using namespace trigger;
0173   // using namespace PhysicsTools;
0174 
0175   test_ev++;
0176   //Get a few collections to apply basic electron ID
0177   //Get electrons
0178   edm::Handle<std::vector<reco::GsfElectron>> electrons;
0179   iEvent.getByToken(electronsToken_, electrons);
0180   //get bdt values
0181   edm::Handle<edm::ValueMap<float>> ele_mva_wp_biased;
0182   iEvent.getByToken(eleBWPToken_, ele_mva_wp_biased);
0183   edm::Handle<edm::ValueMap<float>> ele_mva_wp_unbiased;
0184   iEvent.getByToken(eleUnBWPToken_, ele_mva_wp_unbiased);
0185   //get muons
0186   edm::Handle<std::vector<reco::Muon>> muons;
0187   iEvent.getByToken(muonsToken_, muons);
0188   //Get conversions
0189   edm::Handle<reco::ConversionCollection> conversions;
0190   iEvent.getByToken(conversionsToken_, conversions);
0191   // Get the beam spot
0192   edm::Handle<reco::BeamSpot> theBeamSpot;
0193   iEvent.getByToken(beamSpotToken_, theBeamSpot);
0194   //Get vertices
0195   edm::Handle<std::vector<reco::Vertex>> vertices;
0196   iEvent.getByToken(vtxToken_, vertices);
0197   //continue if there are no vertices
0198   if (vertices->empty())
0199     return false;
0200   edm::Handle<vector<reco::Track>> tracks;
0201   iEvent.getByToken(Tracks_, tracks);
0202   edm::Handle<edm::TriggerResults> trigResults;
0203   iEvent.getByToken(trgresultsToken_, trigResults);
0204   auto const& bField = iSetup.getData(bFieldToken_);
0205   KalmanVertexFitter theKalmanFitter(false);
0206   TransientVertex LLvertex;
0207 
0208   // trigger1=0; trigger2=0; trigger3=0; trigger4=0; trigger5=0; trigger6=0;
0209   nmuons = 0;
0210   nel = 0;
0211   ntracks = 0;
0212 
0213   SelectedMu_index = -1;
0214   SelectedMu_DR = std::numeric_limits<float>::max();
0215   muon_pt.clear();
0216   muon_eta.clear();
0217   muon_phi.clear();
0218   Result = false;
0219   el_pt.clear();
0220   el_eta.clear();
0221   el_phi.clear();
0222   Trk_container.clear();
0223   MuTracks.clear();
0224   ElTracks.clear();
0225   object_container.clear();
0226   object_id.clear();
0227   cleanedTracks.clear();
0228   Epair_ObjectId.clear();
0229   muon_soft.clear();
0230   muon_medium.clear();
0231   muon_tight.clear();
0232   Epair_ObjectIndex.clear();
0233   cleanedObjTracks.clear();
0234   cleanedPairTracks.clear();
0235   Epair_ObjectIndex.clear();
0236   Epair_ObjectId.clear();
0237   Epair_TrkIndex.clear();
0238   //internal stuff
0239   ZvertexTrg = -1 * std::numeric_limits<float>::max();
0240   std::array<float, 5> SelectedTrgObj_PtEtaPhiCharge{{-999, -999, -999, -999, -999}};
0241   vertex_point.SetCoordinates(-1 * std::numeric_limits<float>::max(),
0242                               -1 * std::numeric_limits<float>::max(),
0243                               -1 * std::numeric_limits<float>::max());
0244   for (const reco::Vertex& vtx : *vertices) {
0245     bool isFake = vtx.isFake();
0246     if (isFake)
0247       continue;
0248     vertex_point.SetCoordinates(vtx.x(), vtx.y(), vtx.z());
0249     if (vertex_point.x() != -1 * std::numeric_limits<float>::max())
0250       break;
0251   }
0252   if (vertex_point.x() == -1 * std::numeric_limits<float>::max())
0253     return false;
0254 
0255   beam_x = theBeamSpot->x0();
0256   beam_y = theBeamSpot->y0();
0257   beam_z = theBeamSpot->z0();
0258   if (!hltFired(iEvent, iSetup, HLTPath_))
0259     return false;
0260 
0261   SelectedTrgObj_PtEtaPhiCharge = hltObject(iEvent, iSetup, HLTFilter_);
0262   SelectedMu_DR = std::numeric_limits<float>::max();
0263   MuTracks.clear();
0264   object_container.clear();
0265   object_id.clear();
0266   nmuons = 0;
0267   for (const reco::Muon& mu : *muons) {
0268     if (fabs(mu.eta()) > EtaTrack_Cut)
0269       continue;
0270     //find triggering muon
0271     float deltaRmu = deltaR(mu.eta(), mu.phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]);
0272     if (deltaRmu < MuTrgMatchCone && SelectedMu_DR > deltaRmu) {
0273       SelectedMu_DR = deltaR(mu.eta(), mu.phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]);
0274       ZvertexTrg = mu.vz();
0275     }
0276     //save muons that we want to skim
0277     if (!SkimOnlyElectrons) {
0278       bool tight = false, soft = false;
0279       if (vertices.isValid()) {
0280         tight = isTightMuonCustom(*(&mu), (*vertices)[0]);
0281         soft = muon::isSoftMuon(*(&mu), (*vertices)[0]);
0282       }
0283       const Track* mutrack = mu.bestTrack();
0284       muon_medium.push_back(isMediumMuonCustom(*(&mu)));
0285       muon_tight.push_back(tight);
0286       muon_soft.push_back(soft);
0287       muon_pt.push_back(mu.pt());
0288       muon_eta.push_back(mu.eta());
0289       muon_phi.push_back(mu.phi());
0290       auto muTrack = std::make_shared<reco::Track>(*mutrack);
0291       MuTracks.push_back(muTrack);
0292       object_container.push_back(nmuons);
0293       object_id.push_back(13);
0294       nmuons++;
0295     }
0296   }
0297 
0298   if (SelectedMu_DR == std::numeric_limits<float>::max() && SkipIfNoMuMatch) {
0299     return false;
0300   }
0301 
0302   //Save electrons we want to skim
0303   ElTracks.clear();
0304   if (!SkimOnlyMuons) {
0305     unsigned int count_el = -1;
0306     for (const reco::GsfElectron& el : *electrons) {
0307       count_el++;
0308       bool passConvVeto = !ConversionTools::hasMatchedConversion(*(&el), *conversions, theBeamSpot->position());
0309       if (!passConvVeto)
0310         continue;
0311       reco::GsfTrackRef seed = el.gsfTrack();
0312       if (seed.isNull())
0313         continue;
0314       if ((*ele_mva_wp_biased)[seed] < BiasedWP)
0315         continue;
0316       if ((*ele_mva_wp_unbiased)[seed] < UnbiasedWP)
0317         continue;
0318       const Track* eltrack = el.bestTrack();
0319       if (fabs(eltrack->eta()) > EtaTrack_Cut)
0320         continue;
0321       if (eltrack->pt() < PtEl_Cut)
0322         continue;
0323       auto ElTrack = std::make_shared<reco::Track>(*eltrack);
0324       ElTracks.push_back(ElTrack);
0325       object_container.push_back(nel);
0326       el_pt.push_back(eltrack->pt());
0327       el_eta.push_back(eltrack->eta());
0328       el_phi.push_back(eltrack->phi());
0329       nel++;
0330       object_id.push_back(11);
0331     }
0332   }
0333 
0334   //Save tracks we want to skim: used both as mu or e candidate
0335   cleanedTracks.clear();
0336   trk_index = 0;
0337   for (const reco::Track& trk : *tracks) {
0338     if (!trk.quality(Track::highPurity))
0339       continue;
0340     if (trk.pt() < PtTrack_Cut)
0341       continue;
0342     if (fabs(trk.eta()) > EtaTrack_Cut)
0343       continue;
0344     if (trk.charge() == 0)
0345       continue;
0346     if (trk.normalizedChi2() > MaxChi2Track_Cut || trk.normalizedChi2() < MinChi2Track_Cut)
0347       continue;
0348     if (fabs(trk.dxy()) / trk.dxyError() < TrackSdxy_Cut)
0349       continue;
0350     if (SelectedMu_DR < std::numeric_limits<float>::max()) {
0351       if (fabs(ZvertexTrg - trk.vz()) > TrackMuDz_Cut)
0352         continue;
0353       if (deltaR(trk.eta(), trk.phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0354           TrgExclusionCone)
0355         continue;
0356     }
0357     //assignments
0358     auto cleanTrack = std::make_shared<reco::Track>(trk);
0359     cleanedTracks.push_back(cleanTrack);
0360     Trk_container.push_back(trk_index);
0361     trk_index++;
0362   }
0363 
0364   //create ll combination
0365 
0366   // fit track pairs
0367   cleanedObjTracks.clear();
0368   cleanedPairTracks.clear();
0369   TLorentzVector vel1, vel2;
0370   std::vector<std::shared_ptr<reco::Track>> cleanedObjects;
0371   //add tracks of objects
0372   if (!SkimOnlyElectrons) {
0373     for (auto& vec : MuTracks)
0374       cleanedObjects.push_back(vec);
0375   }
0376 
0377   if (!SkimOnlyMuons) {
0378     for (auto& vec : ElTracks)
0379       cleanedObjects.push_back(vec);
0380   }
0381 
0382   if (cleanedObjects.empty())
0383     return false;
0384 
0385   for (auto& obj : cleanedObjects) {
0386     auto tranobj = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*obj, &bField));
0387     unsigned int iobj = &obj - &cleanedObjects[0];
0388     unsigned int index = object_container.at(iobj);
0389     float massLep = 0.0005;
0390     if (object_id.at(iobj) == 13)
0391       massLep = 0.105;
0392     //posible cut in mu quality
0393     if (object_id.at(iobj) == 13 && QualMu_Cut == 1 && !muon_soft.at(index))
0394       continue;
0395     if (object_id.at(iobj) == 13 && QualMu_Cut == 2 && !muon_medium.at(index))
0396       continue;
0397     if (object_id.at(iobj) == 13 && QualMu_Cut == 3 && !muon_tight.at(index))
0398       continue;
0399     // take the corresponding lepton track for lorentz vector
0400     vel1.SetPtEtaPhiM(obj->pt(), obj->eta(), obj->phi(), massLep);
0401     if (object_id.at(iobj) == 13 && vel1.Pt() < PtMu_Cut)
0402       continue;
0403     for (auto& trk2 : cleanedTracks) {
0404       unsigned int itrk2 = &trk2 - &cleanedTracks[0];
0405       //opposite sign
0406       if (obj->charge() * trk2->charge() == 1)
0407         continue;
0408       if (ObjPtLargerThanTrack && vel1.Pt() < trk2->pt())
0409         continue;
0410       vel2.SetPtEtaPhiM(trk2->pt(), trk2->eta(), trk2->phi(), massLep);
0411       //probe side cuts
0412       if (object_id.at(iobj) == 13 &&
0413           deltaR(obj->eta(), obj->phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0414               MuTrgExclusionCone)
0415         continue;
0416       if (object_id.at(iobj) == 11 &&
0417           deltaR(obj->eta(), obj->phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0418               ElTrgExclusionCone)
0419         continue;
0420       if (SelectedMu_DR < std::numeric_limits<float>::max()) {
0421         if (object_id.at(iobj) == 13 && fabs(ZvertexTrg - obj->vz()) > MuTrgMuDz_Cut)
0422           continue;
0423         if (object_id.at(iobj) == 11 && fabs(ZvertexTrg - obj->vz()) > ElTrgMuDz_Cut)
0424           continue;
0425       }
0426       //additional cuts
0427       float InvMassLepLep = (vel1 + vel2).M();
0428       if (InvMassLepLep > MaxMee_Cut || InvMassLepLep < MinMee_Cut)
0429         continue;
0430       auto trantrk2 = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*trk2, &bField));
0431       std::vector<reco::TransientTrack> tempTracks;
0432       tempTracks.reserve(2);
0433       tempTracks.push_back(*tranobj);
0434       tempTracks.push_back(*trantrk2);
0435       LLvertex = theKalmanFitter.vertex(tempTracks);
0436       if (!LLvertex.isValid())
0437         continue;
0438       if (ChiSquaredProbability(LLvertex.totalChiSquared(), LLvertex.degreesOfFreedom()) < Probee_Cut)
0439         continue;
0440       if (ZvertexTrg > -1 * std::numeric_limits<float>::max() &&
0441           fabs(ZvertexTrg - LLvertex.position().z()) > EpairZvtx_Cut)
0442         continue;
0443       GlobalError err = LLvertex.positionError();
0444       GlobalPoint Dispbeamspot(-1 * ((theBeamSpot->x0() - LLvertex.position().x()) +
0445                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dxdz()),
0446                                -1 * ((theBeamSpot->y0() - LLvertex.position().y()) +
0447                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dydz()),
0448                                0);
0449       math::XYZVector pperp((vel1 + vel2).Px(), (vel1 + vel2).Py(), 0);
0450       math::XYZVector vperp(Dispbeamspot.x(), Dispbeamspot.y(), 0.);
0451       float tempCos = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0452       if (tempCos < Cosee_Cut)
0453         continue;
0454       cleanedObjTracks.push_back(obj);
0455       cleanedPairTracks.push_back(trk2);
0456       Epair_ObjectIndex.push_back(object_container.at(iobj));
0457       Epair_ObjectId.push_back(object_id.at(iobj));
0458       Epair_TrkIndex.push_back(Trk_container.at(itrk2));
0459     }
0460   }
0461 
0462   // B recontrsuction
0463   TLorentzVector vK;
0464   for (unsigned int iobj = 0; iobj < cleanedObjTracks.size(); iobj++) {
0465     auto objtrk = cleanedObjTracks.at(iobj);
0466     auto pairtrk = cleanedPairTracks.at(iobj);
0467     auto tranobj = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*objtrk, &bField));
0468     auto tranpair = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*pairtrk, &bField));
0469     //unsigned int index=Epair_ObjectIndex.at(iobj);
0470     float massLep = 0.0005;
0471     if (Epair_ObjectId.at(iobj) == 13)
0472       massLep = 0.105;
0473     vel1.SetPtEtaPhiM(objtrk->pt(), objtrk->eta(), objtrk->phi(), massLep);
0474     vel2.SetPtEtaPhiM(pairtrk->pt(), pairtrk->eta(), pairtrk->phi(), massLep);
0475     for (auto& trk : cleanedTracks) {
0476       //reject track corresponding to mu or e
0477       if (trk->charge() == objtrk->charge() &&
0478           deltaR(objtrk->eta(), objtrk->phi(), trk->eta(), trk->phi()) < TrkObjExclusionCone)
0479         continue;
0480       if (trk->pt() < PtKTrack_Cut)
0481         continue;
0482       if (fabs(trk->dxy(vertex_point)) / trk->dxyError() < Ksdxy_Cut)
0483         continue;
0484       // skip the track that was used as mu or e
0485       if (Epair_TrkIndex.at(iobj) == &trk - &cleanedTracks[0])
0486         continue;
0487       vK.SetPtEtaPhiM(trk->pt(), trk->eta(), trk->phi(), 0.493);
0488       //final cuts
0489       float InvMass = (vel1 + vel2 + vK).M();
0490       if (InvMass > MaxMB_Cut || InvMass < MinMB_Cut)
0491         continue;
0492       if ((vel1 + vel2 + vK).Pt() < PtB_Cut)
0493         continue;
0494       auto trantrk = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*trk, &bField));
0495       std::vector<reco::TransientTrack> tempTracks;
0496       tempTracks.reserve(3);
0497       tempTracks.push_back(*tranobj);
0498       tempTracks.push_back(*tranpair);
0499       tempTracks.push_back(*trantrk);
0500       LLvertex = theKalmanFitter.vertex(tempTracks);
0501       if (!LLvertex.isValid())
0502         continue;
0503       if (ChiSquaredProbability(LLvertex.totalChiSquared(), LLvertex.degreesOfFreedom()) < ProbeeK_Cut)
0504         continue;
0505       GlobalError err = LLvertex.positionError();
0506       GlobalPoint Dispbeamspot(-1 * ((theBeamSpot->x0() - LLvertex.position().x()) +
0507                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dxdz()),
0508                                -1 * ((theBeamSpot->y0() - LLvertex.position().y()) +
0509                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dydz()),
0510                                0);
0511       math::XYZVector pperp((vel1 + vel2 + vK).Px(), (vel1 + vel2 + vK).Py(), 0);
0512       math::XYZVector vperp(Dispbeamspot.x(), Dispbeamspot.y(), 0.);
0513       float tempCos = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0514       if (tempCos < CoseeK_Cut)
0515         continue;
0516       if (SLxy_Cut > Dispbeamspot.perp() / sqrt(err.rerr(Dispbeamspot)))
0517         continue;
0518       Result = true;
0519       break;
0520     }
0521     if (Result)
0522       break;
0523   }
0524   //decission
0525   return Result;
0526 }
0527 
0528 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0529 void LeptonSkimming::beginStream(edm::StreamID) {}
0530 
0531 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0532 void LeptonSkimming::endStream() {}
0533 
0534 // ------------ method called when starting to processes a run  ------------
0535 /*
0536 void
0537 LeptonSkimming::beginRun(edm::Run const&, edm::EventSetup const&)
0538 { 
0539 }
0540 */
0541 
0542 // ------------ method called when ending the processing of a run  ------------
0543 /*
0544 void
0545 LeptonSkimming::endRun(edm::Run const&, edm::EventSetup const&)
0546 {
0547 }
0548 */
0549 
0550 // ------------ method called when starting to processes a luminosity block  ------------
0551 /*
0552 void
0553 LeptonSkimming::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0554 {
0555 }
0556 */
0557 
0558 // ------------ method called when ending the processing of a luminosity block  ------------
0559 /*
0560 void
0561 LeptonSkimming::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0562 {
0563 }
0564 */
0565 
0566 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0567 void LeptonSkimming::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0568   //The following says we do not know what parameters are allowed so do no validation
0569   // Please change this to state exactly what you do use, even if it is no parameters
0570   edm::ParameterSetDescription desc;
0571   desc.setUnknown();
0572   descriptions.addDefault(desc);
0573 }
0574 //define this as a plug-in
0575 DEFINE_FWK_MODULE(LeptonSkimming);