Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:03:42 

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     for (const reco::GsfElectron& el : *electrons) {
0306       bool passConvVeto = !ConversionTools::hasMatchedConversion(*(&el), *conversions, theBeamSpot->position());
0307       if (!passConvVeto)
0308         continue;
0309       reco::GsfTrackRef seed = el.gsfTrack();
0310       if (seed.isNull())
0311         continue;
0312       if ((*ele_mva_wp_biased)[seed] < BiasedWP)
0313         continue;
0314       if ((*ele_mva_wp_unbiased)[seed] < UnbiasedWP)
0315         continue;
0316       const Track* eltrack = el.bestTrack();
0317       if (fabs(eltrack->eta()) > EtaTrack_Cut)
0318         continue;
0319       if (eltrack->pt() < PtEl_Cut)
0320         continue;
0321       auto ElTrack = std::make_shared<reco::Track>(*eltrack);
0322       ElTracks.push_back(ElTrack);
0323       object_container.push_back(nel);
0324       el_pt.push_back(eltrack->pt());
0325       el_eta.push_back(eltrack->eta());
0326       el_phi.push_back(eltrack->phi());
0327       nel++;
0328       object_id.push_back(11);
0329     }
0330   }
0331 
0332   //Save tracks we want to skim: used both as mu or e candidate
0333   cleanedTracks.clear();
0334   trk_index = 0;
0335   for (const reco::Track& trk : *tracks) {
0336     if (!trk.quality(Track::highPurity))
0337       continue;
0338     if (trk.pt() < PtTrack_Cut)
0339       continue;
0340     if (fabs(trk.eta()) > EtaTrack_Cut)
0341       continue;
0342     if (trk.charge() == 0)
0343       continue;
0344     if (trk.normalizedChi2() > MaxChi2Track_Cut || trk.normalizedChi2() < MinChi2Track_Cut)
0345       continue;
0346     if (fabs(trk.dxy()) / trk.dxyError() < TrackSdxy_Cut)
0347       continue;
0348     if (SelectedMu_DR < std::numeric_limits<float>::max()) {
0349       if (fabs(ZvertexTrg - trk.vz()) > TrackMuDz_Cut)
0350         continue;
0351       if (deltaR(trk.eta(), trk.phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0352           TrgExclusionCone)
0353         continue;
0354     }
0355     //assignments
0356     auto cleanTrack = std::make_shared<reco::Track>(trk);
0357     cleanedTracks.push_back(cleanTrack);
0358     Trk_container.push_back(trk_index);
0359     trk_index++;
0360   }
0361 
0362   //create ll combination
0363 
0364   // fit track pairs
0365   cleanedObjTracks.clear();
0366   cleanedPairTracks.clear();
0367   TLorentzVector vel1, vel2;
0368   std::vector<std::shared_ptr<reco::Track>> cleanedObjects;
0369   //add tracks of objects
0370   if (!SkimOnlyElectrons) {
0371     for (auto& vec : MuTracks)
0372       cleanedObjects.push_back(vec);
0373   }
0374 
0375   if (!SkimOnlyMuons) {
0376     for (auto& vec : ElTracks)
0377       cleanedObjects.push_back(vec);
0378   }
0379 
0380   if (cleanedObjects.empty())
0381     return false;
0382 
0383   for (auto& obj : cleanedObjects) {
0384     auto tranobj = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*obj, &bField));
0385     unsigned int iobj = &obj - &cleanedObjects[0];
0386     unsigned int index = object_container.at(iobj);
0387     float massLep = 0.0005;
0388     if (object_id.at(iobj) == 13)
0389       massLep = 0.105;
0390     //posible cut in mu quality
0391     if (object_id.at(iobj) == 13 && QualMu_Cut == 1 && !muon_soft.at(index))
0392       continue;
0393     if (object_id.at(iobj) == 13 && QualMu_Cut == 2 && !muon_medium.at(index))
0394       continue;
0395     if (object_id.at(iobj) == 13 && QualMu_Cut == 3 && !muon_tight.at(index))
0396       continue;
0397     // take the corresponding lepton track for lorentz vector
0398     vel1.SetPtEtaPhiM(obj->pt(), obj->eta(), obj->phi(), massLep);
0399     if (object_id.at(iobj) == 13 && vel1.Pt() < PtMu_Cut)
0400       continue;
0401     for (auto& trk2 : cleanedTracks) {
0402       unsigned int itrk2 = &trk2 - &cleanedTracks[0];
0403       //opposite sign
0404       if (obj->charge() * trk2->charge() == 1)
0405         continue;
0406       if (ObjPtLargerThanTrack && vel1.Pt() < trk2->pt())
0407         continue;
0408       vel2.SetPtEtaPhiM(trk2->pt(), trk2->eta(), trk2->phi(), massLep);
0409       //probe side cuts
0410       if (object_id.at(iobj) == 13 &&
0411           deltaR(obj->eta(), obj->phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0412               MuTrgExclusionCone)
0413         continue;
0414       if (object_id.at(iobj) == 11 &&
0415           deltaR(obj->eta(), obj->phi(), SelectedTrgObj_PtEtaPhiCharge[1], SelectedTrgObj_PtEtaPhiCharge[2]) <
0416               ElTrgExclusionCone)
0417         continue;
0418       if (SelectedMu_DR < std::numeric_limits<float>::max()) {
0419         if (object_id.at(iobj) == 13 && fabs(ZvertexTrg - obj->vz()) > MuTrgMuDz_Cut)
0420           continue;
0421         if (object_id.at(iobj) == 11 && fabs(ZvertexTrg - obj->vz()) > ElTrgMuDz_Cut)
0422           continue;
0423       }
0424       //additional cuts
0425       float InvMassLepLep = (vel1 + vel2).M();
0426       if (InvMassLepLep > MaxMee_Cut || InvMassLepLep < MinMee_Cut)
0427         continue;
0428       auto trantrk2 = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*trk2, &bField));
0429       std::vector<reco::TransientTrack> tempTracks;
0430       tempTracks.reserve(2);
0431       tempTracks.push_back(*tranobj);
0432       tempTracks.push_back(*trantrk2);
0433       LLvertex = theKalmanFitter.vertex(tempTracks);
0434       if (!LLvertex.isValid())
0435         continue;
0436       if (ChiSquaredProbability(LLvertex.totalChiSquared(), LLvertex.degreesOfFreedom()) < Probee_Cut)
0437         continue;
0438       if (ZvertexTrg > -1 * std::numeric_limits<float>::max() &&
0439           fabs(ZvertexTrg - LLvertex.position().z()) > EpairZvtx_Cut)
0440         continue;
0441       GlobalError err = LLvertex.positionError();
0442       GlobalPoint Dispbeamspot(-1 * ((theBeamSpot->x0() - LLvertex.position().x()) +
0443                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dxdz()),
0444                                -1 * ((theBeamSpot->y0() - LLvertex.position().y()) +
0445                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dydz()),
0446                                0);
0447       math::XYZVector pperp((vel1 + vel2).Px(), (vel1 + vel2).Py(), 0);
0448       math::XYZVector vperp(Dispbeamspot.x(), Dispbeamspot.y(), 0.);
0449       float tempCos = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0450       if (tempCos < Cosee_Cut)
0451         continue;
0452       cleanedObjTracks.push_back(obj);
0453       cleanedPairTracks.push_back(trk2);
0454       Epair_ObjectIndex.push_back(object_container.at(iobj));
0455       Epair_ObjectId.push_back(object_id.at(iobj));
0456       Epair_TrkIndex.push_back(Trk_container.at(itrk2));
0457     }
0458   }
0459 
0460   // B recontrsuction
0461   TLorentzVector vK;
0462   for (unsigned int iobj = 0; iobj < cleanedObjTracks.size(); iobj++) {
0463     auto objtrk = cleanedObjTracks.at(iobj);
0464     auto pairtrk = cleanedPairTracks.at(iobj);
0465     auto tranobj = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*objtrk, &bField));
0466     auto tranpair = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*pairtrk, &bField));
0467     //unsigned int index=Epair_ObjectIndex.at(iobj);
0468     float massLep = 0.0005;
0469     if (Epair_ObjectId.at(iobj) == 13)
0470       massLep = 0.105;
0471     vel1.SetPtEtaPhiM(objtrk->pt(), objtrk->eta(), objtrk->phi(), massLep);
0472     vel2.SetPtEtaPhiM(pairtrk->pt(), pairtrk->eta(), pairtrk->phi(), massLep);
0473     for (auto& trk : cleanedTracks) {
0474       //reject track corresponding to mu or e
0475       if (trk->charge() == objtrk->charge() &&
0476           deltaR(objtrk->eta(), objtrk->phi(), trk->eta(), trk->phi()) < TrkObjExclusionCone)
0477         continue;
0478       if (trk->pt() < PtKTrack_Cut)
0479         continue;
0480       if (fabs(trk->dxy(vertex_point)) / trk->dxyError() < Ksdxy_Cut)
0481         continue;
0482       // skip the track that was used as mu or e
0483       if (Epair_TrkIndex.at(iobj) == &trk - &cleanedTracks[0])
0484         continue;
0485       vK.SetPtEtaPhiM(trk->pt(), trk->eta(), trk->phi(), 0.493);
0486       //final cuts
0487       float InvMass = (vel1 + vel2 + vK).M();
0488       if (InvMass > MaxMB_Cut || InvMass < MinMB_Cut)
0489         continue;
0490       if ((vel1 + vel2 + vK).Pt() < PtB_Cut)
0491         continue;
0492       auto trantrk = std::make_shared<reco::TransientTrack>(reco::TransientTrack(*trk, &bField));
0493       std::vector<reco::TransientTrack> tempTracks;
0494       tempTracks.reserve(3);
0495       tempTracks.push_back(*tranobj);
0496       tempTracks.push_back(*tranpair);
0497       tempTracks.push_back(*trantrk);
0498       LLvertex = theKalmanFitter.vertex(tempTracks);
0499       if (!LLvertex.isValid())
0500         continue;
0501       if (ChiSquaredProbability(LLvertex.totalChiSquared(), LLvertex.degreesOfFreedom()) < ProbeeK_Cut)
0502         continue;
0503       GlobalError err = LLvertex.positionError();
0504       GlobalPoint Dispbeamspot(-1 * ((theBeamSpot->x0() - LLvertex.position().x()) +
0505                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dxdz()),
0506                                -1 * ((theBeamSpot->y0() - LLvertex.position().y()) +
0507                                      (LLvertex.position().z() - theBeamSpot->z0()) * theBeamSpot->dydz()),
0508                                0);
0509       math::XYZVector pperp((vel1 + vel2 + vK).Px(), (vel1 + vel2 + vK).Py(), 0);
0510       math::XYZVector vperp(Dispbeamspot.x(), Dispbeamspot.y(), 0.);
0511       float tempCos = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0512       if (tempCos < CoseeK_Cut)
0513         continue;
0514       if (SLxy_Cut > Dispbeamspot.perp() / sqrt(err.rerr(Dispbeamspot)))
0515         continue;
0516       Result = true;
0517       break;
0518     }
0519     if (Result)
0520       break;
0521   }
0522   //decission
0523   return Result;
0524 }
0525 
0526 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0527 void LeptonSkimming::beginStream(edm::StreamID) {}
0528 
0529 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0530 void LeptonSkimming::endStream() {}
0531 
0532 // ------------ method called when starting to processes a run  ------------
0533 /*
0534 void
0535 LeptonSkimming::beginRun(edm::Run const&, edm::EventSetup const&)
0536 { 
0537 }
0538 */
0539 
0540 // ------------ method called when ending the processing of a run  ------------
0541 /*
0542 void
0543 LeptonSkimming::endRun(edm::Run const&, edm::EventSetup const&)
0544 {
0545 }
0546 */
0547 
0548 // ------------ method called when starting to processes a luminosity block  ------------
0549 /*
0550 void
0551 LeptonSkimming::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0552 {
0553 }
0554 */
0555 
0556 // ------------ method called when ending the processing of a luminosity block  ------------
0557 /*
0558 void
0559 LeptonSkimming::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0560 {
0561 }
0562 */
0563 
0564 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0565 void LeptonSkimming::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0566   //The following says we do not know what parameters are allowed so do no validation
0567   // Please change this to state exactly what you do use, even if it is no parameters
0568   edm::ParameterSetDescription desc;
0569   desc.setUnknown();
0570   descriptions.addDefault(desc);
0571 }
0572 //define this as a plug-in
0573 DEFINE_FWK_MODULE(LeptonSkimming);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!