Reference Manual

Go to the documentation of this file.
00001 //
00002 // $Id: TtDilepEvtSolutionMaker.cc,v 1.26 2010/04/30 12:52:20 dammann Exp $
00003 //
00004 
00005 #include "DataFormats/Math/interface/deltaR.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 
00008 #include "TopQuarkAnalysis/TopEventSelection/interface/TtDilepLRSignalSelObservables.h"
00009 
00010 #include "AnalysisDataFormats/TopObjects/interface/TtDilepEvtSolution.h"
00011 #include "TopQuarkAnalysis/TopEventProducers/interface/TtDilepEvtSolutionMaker.h"
00012 
00013 #include <memory>
00014 #include <vector>
00015 
00016 
00017 TtDilepEvtSolutionMaker::TtDilepEvtSolutionMaker(const edm::ParameterSet & iConfig) 
00018 {
00019   // configurables
00020   electronSource_ = iConfig.getParameter<edm::InputTag>("electronSource");
00021   muonSource_     = iConfig.getParameter<edm::InputTag>("muonSource");
00022   tauSource_      = iConfig.getParameter<edm::InputTag>("tauSource");
00023   metSource_      = iConfig.getParameter<edm::InputTag>("metSource");
00024   jetSource_      = iConfig.getParameter<edm::InputTag>("jetSource");
00025   jetCorrScheme_  = iConfig.getParameter<int>          ("jetCorrectionScheme");
00026   evtSource_      = iConfig.getParameter<edm::InputTag>("evtSource");
00027   nrCombJets_     = iConfig.getParameter<unsigned int> ("nrCombJets");
00028   matchToGenEvt_  = iConfig.getParameter<bool>         ("matchToGenEvt");
00029   calcTopMass_    = iConfig.getParameter<bool>         ("calcTopMass"); 
00030   useMCforBest_   = iConfig.getParameter<bool>         ("bestSolFromMC");
00031   eeChannel_      = iConfig.getParameter<bool>         ("eeChannel"); 
00032   emuChannel_     = iConfig.getParameter<bool>         ("emuChannel");
00033   mumuChannel_    = iConfig.getParameter<bool>         ("mumuChannel");
00034   mutauChannel_   = iConfig.getParameter<bool>         ("mutauChannel");
00035   etauChannel_    = iConfig.getParameter<bool>         ("etauChannel");
00036   tautauChannel_  = iConfig.getParameter<bool>         ("tautauChannel");
00037   tmassbegin_     = iConfig.getParameter<double>       ("tmassbegin");
00038   tmassend_       = iConfig.getParameter<double>       ("tmassend");
00039   tmassstep_      = iConfig.getParameter<double>       ("tmassstep");
00040   nupars_         = iConfig.getParameter<std::vector<double> >("neutrino_parameters");
00041     
00042   // define what will be produced
00043   produces<std::vector<TtDilepEvtSolution> >();
00044   
00045   myLRSignalSelObservables = new TtDilepLRSignalSelObservables();
00046   myLRSignalSelObservables->jetSource(jetSource_);
00047 }
00048 
00049 TtDilepEvtSolutionMaker::~TtDilepEvtSolutionMaker() 
00050 {
00051 }
00052 
00053 void TtDilepEvtSolutionMaker::beginJob()
00054 {
00055   solver = new TtFullLepKinSolver(tmassbegin_, tmassend_, tmassstep_, nupars_);
00056 }
00057 
00058 void TtDilepEvtSolutionMaker::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) 
00059 {
00060   edm::Handle<std::vector<pat::Tau> > taus;
00061   iEvent.getByLabel(tauSource_, taus);
00062   edm::Handle<std::vector<pat::Muon> > muons;
00063   iEvent.getByLabel(muonSource_, muons);
00064   edm::Handle<std::vector<pat::Electron> > electrons;
00065   iEvent.getByLabel(electronSource_, electrons);
00066   edm::Handle<std::vector<pat::MET> > mets;
00067   iEvent.getByLabel(metSource_, mets);
00068   edm::Handle<std::vector<pat::Jet> > jets;
00069   iEvent.getByLabel(jetSource_, jets);
00070   
00071   int selMuonp = -1, selMuonm = -1;
00072   int selElectronp = -1, selElectronm = -1;
00073   int selTaup = -1, selTaum = -1;
00074   bool leptonFound = false;
00075   bool mumu = false;
00076   bool emu = false;
00077   bool ee = false;
00078   bool etau = false;
00079   bool mutau = false;
00080   bool tautau = false;
00081   bool leptonFoundEE = false;
00082   bool leptonFoundMM = false;
00083   bool leptonFoundTT = false;
00084   bool leptonFoundEpMm = false;
00085   bool leptonFoundEmMp = false;
00086   bool leptonFoundEpTm = false;
00087   bool leptonFoundEmTp = false;
00088   bool leptonFoundMpTm = false;
00089   bool leptonFoundMmTp = false;
00090   bool jetsFound = false;
00091   bool METFound = false;
00092   std::vector<int>  JetVetoByTaus;
00093   
00094   //select MET (TopMET vector is sorted on ET)
00095   if(mets->size()>=1) { METFound = true; }
00096   
00097   // If we have electrons and muons available, 
00098   // build a solutions with electrons and muons.
00099   if (muons->size() + electrons->size() >=2) {
00100     // select leptons
00101     if (electrons->size() == 0) mumu = true;
00102     else if (muons->size() == 0) ee = true;
00103     else if (electrons->size() == 1) {
00104       if (muons->size() == 1) emu = true;
00105       else if (PTComp(&(*electrons)[0], &(*muons)[1])) emu = true;
00106       else  mumu = true;
00107     }
00108     else if (electrons->size() > 1) {
00109       if (PTComp(&(*electrons)[1], &(*muons)[0])) ee = true;
00110       else if (muons->size() == 1) emu = true;
00111       else if (PTComp(&(*electrons)[0], &(*muons)[1])) emu = true;
00112       else mumu = true;
00113     }
00114     if (ee) {
00115       if (LepDiffCharge(&(*electrons)[0], &(*electrons)[1])) {
00116         leptonFound = true;
00117         leptonFoundEE = true;
00118         if (HasPositiveCharge(&(*electrons)[0])) {
00119           selElectronp = 0;
00120           selElectronm = 1;
00121         } else {
00122           selElectronp = 1;
00123           selElectronm = 0;
00124         }
00125       }
00126     }
00127     else if (emu) {
00128       if (LepDiffCharge(&(*electrons)[0], &(*muons)[0])) {
00129         leptonFound = true;
00130         if (HasPositiveCharge(&(*electrons)[0])) {
00131           leptonFoundEpMm = true;
00132           selElectronp = 0;
00133           selMuonm = 0;
00134         } else {
00135           leptonFoundEmMp = true;
00136           selMuonp = 0;
00137           selElectronm = 0;
00138         }
00139       }
00140     }
00141     else if (mumu) {
00142       if (LepDiffCharge(&(*muons)[0], &(*muons)[1])) {
00143         leptonFound = true;
00144         leptonFoundMM = true;
00145         if (HasPositiveCharge(&(*muons)[0])) {
00146           selMuonp = 0;
00147           selMuonm = 1;
00148         } else {
00149           selMuonp = 1;
00150           selMuonm = 0;
00151         }
00152       }
00153     }
00154     //select Jets (TopJet vector is sorted on ET)
00155     if(jets->size()>=2) { jetsFound = true; }
00156   }
00157   // If a tau is needed to have two leptons, then only consider the taus.
00158   // This is the minimal modification of the dilept selection that includes taus,
00159   // since we are considering taus only when no other solution exist.
00160   else if(muons->size() + electrons->size()==1 && taus->size()>0) {
00161     // select leptons
00162     if(muons->size()==1) {
00163       mutau = true;
00164       // depending on the muon charge, set the right muon index and specify channel
00165       int expectedCharge = - muons->begin()->charge();
00166       int* tauIdx = NULL;
00167       if (expectedCharge<0) {
00168         selMuonp = 0;
00169         tauIdx = &selTaum;
00170         leptonFoundMpTm = true;
00171       } else {
00172         selMuonm = 0;
00173         tauIdx = &selTaup;
00174         leptonFoundMmTp = true;
00175       }
00176       // loop over the vector of taus to find the ones
00177       // that have the charge opposite to the muon one, and do not match in eta-phi
00178       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
00179       for(std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau ) {
00180         if(tau->charge()*expectedCharge>=0 && DeltaR<pat::Particle>()(*tau,*(muons->begin()))>0.1) { 
00181           *tauIdx = tau-taus->begin(); 
00182           leptonFound = true;
00183           subset1.push_back(tau);
00184         }
00185       }
00186       // if there are more than one tau with ecalIsol==0, take the smallest E/P
00187       float iso = 999.;
00188       for(std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin(); tau < subset1.end(); ++tau) {
00189         if((*tau)->isCaloTau() && (*tau)->isolationTracksPtSum()<iso) {
00190           *tauIdx = *tau - taus->begin();
00191           iso = (*tau)->isolationTracksPtSum();
00192         }
00193         if((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum()<iso) {
00194           *tauIdx = *tau - taus->begin();
00195           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
00196         }
00197       }
00198       
00199       // check that one combination has been found
00200       if(!leptonFound) { leptonFoundMpTm = false; leptonFoundMmTp = false; } 
00201       // discard the jet that matches the tau (if one) 
00202       if(leptonFound) {
00203         for(std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet<jets->end(); ++jet) {
00204           if(DeltaR<pat::Particle, pat::Jet>()(*(taus->begin()+*tauIdx),*jet)<0.1) {
00205             JetVetoByTaus.push_back(jet-jets->begin());
00206           }
00207         }
00208       }
00209     }
00210     else {
00211       etau = true;
00212       // depending on the electron charge, set the right electron index and specify channel
00213       int expectedCharge = - electrons->begin()->charge();
00214       int* tauIdx = NULL;
00215       if (expectedCharge<0) {
00216         selElectronp = 0;
00217         tauIdx = &selTaum;
00218         leptonFoundEpTm = true;
00219       } else {
00220         selElectronm = 0;
00221         tauIdx = &selTaup;
00222         leptonFoundEmTp = true;
00223       }
00224       // loop over the vector of taus to find the ones
00225       // that have the charge opposite to the muon one, and do not match in eta-phi
00226       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
00227       for(std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau ) {
00228         if(tau->charge()*expectedCharge>=0 && DeltaR<pat::Particle>()(*tau,*(electrons->begin()))>0.1) { 
00229           *tauIdx = tau-taus->begin(); 
00230           leptonFound = true; 
00231           subset1.push_back(tau);
00232         }
00233       }
00234       // if there are more than one tau with ecalIsol==0, take the smallest E/P
00235       float iso = 999.;
00236       for(std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin(); tau < subset1.end(); ++tau) {
00237         if((*tau)->isCaloTau() && (*tau)->isolationTracksPtSum()<iso) {
00238           *tauIdx = *tau - taus->begin();
00239           iso = (*tau)->isolationTracksPtSum();
00240         }
00241         if((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum()<iso) {
00242           *tauIdx = *tau - taus->begin();
00243           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
00244         }
00245       }
00246 
00247       // check that one combination has been found
00248       if(!leptonFound) { leptonFoundEpTm = false; leptonFoundEmTp = false; } 
00249       // discard the jet that matches the tau (if one) 
00250       if(leptonFound) {
00251         for(std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet<jets->end(); ++jet) {
00252           if(DeltaR<pat::Particle, pat::Jet>()(*(taus->begin()+*tauIdx),*jet)<0.1) {
00253             JetVetoByTaus.push_back(jet-jets->begin());
00254           }
00255         }
00256       }
00257     }
00258     // select Jets (TopJet vector is sorted on ET)
00259     jetsFound = ((jets->size()-JetVetoByTaus.size())>=2);
00260   } else if(taus->size()>1) {
00261     tautau = true;
00262     if(LepDiffCharge(&(*taus)[0],&(*taus)[1])) {
00263       leptonFound = true;
00264       leptonFoundTT = true;
00265       if(HasPositiveCharge(&(*taus)[0])) {
00266         selTaup = 0;
00267         selTaum = 1;
00268       }
00269       else {
00270         selTaup = 1;
00271         selTaum = 0;
00272       }
00273     }
00274     for(std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet<jets->end(); ++jet) {
00275       if(DeltaR<pat::Particle, pat::Jet>()((*taus)[0],*jet)<0.1 || DeltaR<pat::Particle, pat::Jet>()((*taus)[1],*jet)<0.1) {
00276         JetVetoByTaus.push_back(jet-jets->begin());
00277       }
00278     }
00279     // select Jets (TopJet vector is sorted on ET)
00280     jetsFound = ((jets->size()-JetVetoByTaus.size())>=2);
00281   }
00282  
00283   // Check that the above work makes sense
00284   if(int(ee)+int(emu)+int(mumu)+int(etau)+int(mutau)+int(tautau)>1) 
00285     std::cout << "[TtDilepEvtSolutionMaker]: "
00286               << "Lepton selection criteria uncorrectly defined" << std::endl;
00287   
00288   bool correctLepton = (leptonFoundEE && eeChannel_)                          ||
00289                        ((leptonFoundEmMp || leptonFoundEpMm) && emuChannel_)  ||
00290                        (leptonFoundMM && mumuChannel_)                        ||
00291                        ((leptonFoundMmTp || leptonFoundMpTm) && mutauChannel_)||
00292                        ((leptonFoundEmTp || leptonFoundEpTm) && etauChannel_) ||
00293                        (leptonFoundTT && tautauChannel_)                        ;
00294                        
00295   std::vector<TtDilepEvtSolution> * evtsols = new std::vector<TtDilepEvtSolution>();
00296   if(correctLepton && METFound && jetsFound) {
00297     // protect against reading beyond array boundaries while discounting vetoed jets
00298     unsigned int nrCombJets = 0; 
00299     unsigned int numberOfJets = 0;
00300     for(; nrCombJets<jets->size() && numberOfJets<nrCombJets_; ++nrCombJets) {
00301       if(find(JetVetoByTaus.begin(),JetVetoByTaus.end(),int(nrCombJets))==JetVetoByTaus.end()) ++numberOfJets;
00302     }
00303     // consider all permutations
00304     for (unsigned int ib = 0; ib < nrCombJets; ib++) {
00305       // skipped jet vetoed during components-flagging.
00306       if(find(JetVetoByTaus.begin(),JetVetoByTaus.end(),int(ib))!=JetVetoByTaus.end())continue;
00307       // second loop of the permutations
00308       for (unsigned int ibbar = 0; ibbar < nrCombJets; ibbar++) {
00309         // avoid the diagonal: b and bbar must be distinct jets
00310         if(ib==ibbar) continue;
00311         // skipped jet vetoed during components-flagging.
00312         if(find(JetVetoByTaus.begin(),JetVetoByTaus.end(),int(ibbar))!=JetVetoByTaus.end())continue;
00313         // Build and save a solution
00314         TtDilepEvtSolution asol;
00315         asol.setJetCorrectionScheme(jetCorrScheme_);
00316         double xconstraint = 0, yconstraint = 0;
00317         // Set e+ in the event
00318         if (leptonFoundEE || leptonFoundEpMm || leptonFoundEpTm) {
00319           asol.setElectronp(electrons, selElectronp);
00320           xconstraint += (*electrons)[selElectronp].px();
00321           yconstraint += (*electrons)[selElectronp].py();
00322         }
00323         // Set e- in the event
00324         if (leptonFoundEE || leptonFoundEmMp || leptonFoundEmTp) {
00325           asol.setElectronm(electrons, selElectronm);
00326           xconstraint += (*electrons)[selElectronm].px();
00327           yconstraint += (*electrons)[selElectronm].py();
00328         }
00329         // Set mu+ in the event
00330         if (leptonFoundMM || leptonFoundEmMp || leptonFoundMpTm) {
00331           asol.setMuonp(muons, selMuonp);
00332           xconstraint += (*muons)[selMuonp].px();
00333           yconstraint += (*muons)[selMuonp].py();
00334         }
00335         // Set mu- in the event
00336         if (leptonFoundMM || leptonFoundEpMm || leptonFoundMmTp) {
00337           asol.setMuonm(muons, selMuonm);
00338           xconstraint += (*muons)[selMuonm].px();
00339           yconstraint += (*muons)[selMuonm].py();
00340         }
00341         // Set tau- in the event
00342         if (leptonFoundEpTm || leptonFoundMpTm || leptonFoundTT) {
00343           asol.setTaum(taus, selTaum);
00344           xconstraint += (*taus)[selTaum].px();
00345           yconstraint += (*taus)[selTaum].py();
00346         }
00347         // Set tau+ in the event
00348         if (leptonFoundEmTp || leptonFoundMmTp || leptonFoundTT) {
00349           asol.setTaup(taus, selTaup);
00350           xconstraint += (*taus)[selTaup].px();
00351           yconstraint += (*taus)[selTaup].py();
00352         }
00353         // Set Jets/MET in the event
00354         asol.setB(jets, ib); 
00355         asol.setBbar(jets, ibbar);
00356         asol.setMET(mets, 0);
00357         xconstraint += (*jets)[ib].px() + (*jets)[ibbar].px() + (*mets)[0].px();
00358         yconstraint += (*jets)[ib].py() + (*jets)[ibbar].py() + (*mets)[0].py();
00359         // if asked for, match the event solutions to the gen Event
00360         if(matchToGenEvt_){
00361           edm::Handle<TtGenEvent> genEvt;
00362           iEvent.getByLabel (evtSource_,genEvt);
00363           asol.setGenEvt(genEvt);
00364         } 
00365         // If asked, use the kin fitter to compute the top mass
00366         if (calcTopMass_) {
00367           solver->SetConstraints(xconstraint, yconstraint);
00368           solver->useWeightFromMC(useMCforBest_);
00369           asol = solver->addKinSolInfo(&asol);
00370         }
00371 
00372      // these lines calculate the observables to be used in the TtDilepSignalSelection LR
00373       (*myLRSignalSelObservables)(asol, iEvent);
00374 
00375         evtsols->push_back(asol);
00376       }
00377     } 
00378     // flag the best solution (MC matching)
00379     if(matchToGenEvt_){
00380       double bestSolDR = 9999.;
00381       int bestSol = -1;
00382       double dR = 0.;
00383       for(size_t s=0; s<evtsols->size(); s++) {
00384         dR = (*evtsols)[s].getJetResidual();
00385         if(dR<bestSolDR) { bestSolDR = dR; bestSol = s; }
00386       }
00387       if(bestSol!=-1) (*evtsols)[bestSol].setBestSol(true);
00388     }
00389     // put the result in the event
00390     std::auto_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
00391     iEvent.put(pOut);
00392   } else {
00393     // no solution: put a dummy solution in the event
00394     TtDilepEvtSolution asol;
00395     evtsols->push_back(asol);
00396     std::auto_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
00397     iEvent.put(pOut);
00398   }
00399 }
00400
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/TopQuarkAnalysis/TopEventProducers/src/TtDilepEvtSolutionMaker.cc
