CMS 3D CMS Logo

HLTTauMCProducer.cc
Go to the documentation of this file.
2 
3 using namespace edm;
4 using namespace std;
5 using namespace reco;
6 
8  : MC_{consumes<GenParticleCollection>(mc.getUntrackedParameter<edm::InputTag>("GenParticles"))},
9  MCMET_{consumes<GenMETCollection>(mc.getUntrackedParameter<edm::InputTag>("GenMET"))},
10  ptMinMCTau_{mc.getUntrackedParameter<double>("ptMinTau", 5.)},
11  ptMinMCElectron_{mc.getUntrackedParameter<double>("ptMinElectron", 5.)},
12  ptMinMCMuon_{mc.getUntrackedParameter<double>("ptMinMuon", 2.)},
13  m_PDG_{mc.getUntrackedParameter<std::vector<int>>("BosonID")},
14  etaMin_{mc.getUntrackedParameter<double>("EtaMin", -2.5)},
15  etaMax_{mc.getUntrackedParameter<double>("EtaMax", 2.5)},
16  phiMin_{mc.getUntrackedParameter<double>("PhiMin", -3.15)},
17  phiMax_{mc.getUntrackedParameter<double>("PhiMax", 3.15)} {
18  // One Parameter Set per Collection
19 
20  produces<LorentzVectorCollection>("LeptonicTauLeptons");
21  produces<LorentzVectorCollection>("LeptonicTauElectrons");
22  produces<LorentzVectorCollection>("LeptonicTauMuons");
23  produces<LorentzVectorCollection>("HadronicTauOneProng");
24  produces<LorentzVectorCollection>("HadronicTauThreeProng");
25  produces<LorentzVectorCollection>("HadronicTauOneAndThreeProng");
26  produces<LorentzVectorCollection>("TauOther");
27  produces<LorentzVectorCollection>("Neutrina");
28  produces<LorentzVectorCollection>("MET");
29  produces<std::vector<int>>("Mothers");
30 }
31 
33  // All the code from HLTTauMCInfo is here :-)
34 
35  unique_ptr<LorentzVectorCollection> product_Electrons(new LorentzVectorCollection);
36  unique_ptr<LorentzVectorCollection> product_Muons(new LorentzVectorCollection);
37  unique_ptr<LorentzVectorCollection> product_Leptons(new LorentzVectorCollection);
38  unique_ptr<LorentzVectorCollection> product_OneProng(new LorentzVectorCollection);
39  unique_ptr<LorentzVectorCollection> product_ThreeProng(new LorentzVectorCollection);
40  unique_ptr<LorentzVectorCollection> product_OneAndThreeProng(new LorentzVectorCollection);
41  unique_ptr<LorentzVectorCollection> product_Other(new LorentzVectorCollection);
42  unique_ptr<LorentzVectorCollection> product_Neutrina(new LorentzVectorCollection);
43  unique_ptr<LorentzVectorCollection> product_MET(new LorentzVectorCollection);
44  unique_ptr<std::vector<int>> product_Mothers(new std::vector<int>);
45 
47  iEvent.getByToken(MC_, genParticles);
48 
49  if (!genParticles.isValid())
50  return;
51 
52  // Look for MET
54  iEvent.getByToken(MCMET_, genMet);
55  LorentzVector MET(0., 0., 0., 0.);
56  if (genMet.isValid()) {
57  MET = LorentzVector(genMet->front().px(), genMet->front().py(), 0, genMet->front().pt());
58  }
59  product_MET->push_back(MET);
60 
61  // Look for primary bosons
62  // It is not guaranteed that primary bosons are stored in event history.
63  // Is it really needed when check if taus from the boson is removed?
64  // Kept for backward compatibility
65  for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p) {
66  // Check the PDG ID
67  bool pdg_ok = false;
68  for (size_t pi = 0; pi < m_PDG_.size(); ++pi) {
69  if (abs((*p).pdgId()) == m_PDG_[pi] && ((*p).isHardProcess() || (*p).status() == 3)) {
70  pdg_ok = true;
71  // cout<<" Bsoson particles: "<< (*p).pdgId()<< " " <<(*p).status() << "
72  // "<< pdg_ok<<endl;
73  break;
74  }
75  }
76 
77  // Check if the boson is one of interest and if there is a valid vertex
78  if (pdg_ok) {
79  product_Mothers->push_back((*p).pdgId());
80 
81  TLorentzVector Boson((*p).px(), (*p).py(), (*p).pz(), (*p).energy());
82  }
83  } // End of search for the bosons
84 
85  // Look for taus
86  GenParticleRefVector allTaus;
87  unsigned index = 0;
88  for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p, ++index) {
89  const GenParticle &genP = *p;
90  // accept only isPromptDecayed() particles
91  if (!genP.isPromptDecayed())
92  continue;
93  // check if it is tau, i.e. if |pdgId|=15
94  if (std::abs(genP.pdgId()) == 15) {
96  // check if it is the last tau in decay/radiation chain
97  GenParticleRefVector daugTaus;
98  getGenDecayProducts(genRef, daugTaus, 0, 15);
99  if (daugTaus.empty())
100  allTaus.push_back(genRef);
101  }
102  }
103 
104  // Find stable tau decay products and build visible taus
105  for (GenParticleRefVector::const_iterator t = allTaus.begin(); t != allTaus.end(); ++t) {
106  // look for all stable (status=1) decay products
107  GenParticleRefVector decayProducts;
108  getGenDecayProducts(*t, decayProducts, 1);
109 
110  // build visible taus and recognize decay mode
111  if (!decayProducts.empty()) {
112  LorentzVector Visible_Taus(0., 0., 0., 0.);
113  LorentzVector TauDecayProduct(0., 0., 0., 0.);
114  LorentzVector Neutrino(0., 0., 0., 0.);
115 
116  int numElectrons = 0;
117  int numMuons = 0;
118  int numChargedPions = 0;
119  int numNeutralPions = 0;
120  int numPhotons = 0;
121  int numNeutrinos = 0;
122  int numOtherParticles = 0;
123 
124  for (GenParticleRefVector::const_iterator pit = decayProducts.begin(); pit != decayProducts.end(); ++pit) {
125  int pdg_id = abs((*pit)->pdgId());
126  if (pdg_id == 11)
127  numElectrons++;
128  else if (pdg_id == 13)
129  numMuons++;
130  else if (pdg_id == 211 || pdg_id == 321)
131  numChargedPions++; // Count both pi+ and K+
132  else if (pdg_id == 111 || pdg_id == 130 || pdg_id == 310)
133  numNeutralPions++; // Count both pi0 and K0_L/S
134  else if (pdg_id == 12 || pdg_id == 14 || pdg_id == 16) {
135  numNeutrinos++;
136  if (pdg_id == 16) {
137  Neutrino.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
138  }
139  } else if (pdg_id == 22)
140  numPhotons++;
141  else {
142  numOtherParticles++;
143  }
144 
145  if (pdg_id != 12 && pdg_id != 14 && pdg_id != 16) {
146  TauDecayProduct.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
147  Visible_Taus += TauDecayProduct;
148  }
149  // cout<< "This has to be the same: " << (*pit)->pdgId()
150  //<< " "<< (*pit)->status()<< " mother: "<< (*pit)->mother()->pdgId() <<
151  // endl;
152  }
153 
154  int tauDecayMode = kOther;
155 
156  if (numOtherParticles == 0) {
157  if (numElectrons == 1) {
158  //--- tau decays into electrons
159  tauDecayMode = kElectron;
160  } else if (numMuons == 1) {
161  //--- tau decays into muons
162  tauDecayMode = kMuon;
163  } else {
164  //--- hadronic tau decays
165  switch (numChargedPions) {
166  case 1:
167  if (numNeutralPions != 0) {
168  tauDecayMode = kOther;
169  break;
170  }
171  switch (numPhotons) {
172  case 0:
173  tauDecayMode = kOneProng0pi0;
174  break;
175  case 2:
176  tauDecayMode = kOneProng1pi0;
177  break;
178  case 4:
179  tauDecayMode = kOneProng2pi0;
180  break;
181  default:
182  tauDecayMode = kOther;
183  break;
184  }
185  break;
186  case 3:
187  if (numNeutralPions != 0) {
188  tauDecayMode = kOther;
189  break;
190  }
191  switch (numPhotons) {
192  case 0:
193  tauDecayMode = kThreeProng0pi0;
194  break;
195  case 2:
196  tauDecayMode = kThreeProng1pi0;
197  break;
198  default:
199  tauDecayMode = kOther;
200  break;
201  }
202  break;
203  }
204  }
205  }
206 
207  // cout<< "So we have a: " << tauDecayMode <<endl;
208  if (tauDecayMode == kElectron) {
209  if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
210  Visible_Taus.phi() < phiMax_) &&
211  (Visible_Taus.pt() > ptMinMCElectron_)) {
212  product_Electrons->push_back(Visible_Taus);
213  product_Leptons->push_back(Visible_Taus);
214  }
215  } else if (tauDecayMode == kMuon) {
216  if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
217  Visible_Taus.phi() < phiMax_) &&
218  (Visible_Taus.pt() > ptMinMCMuon_)) {
219  product_Muons->push_back(Visible_Taus);
220  product_Leptons->push_back(Visible_Taus);
221  }
222  } else if (tauDecayMode == kOneProng0pi0 || tauDecayMode == kOneProng1pi0 || tauDecayMode == kOneProng2pi0) {
223  if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
224  Visible_Taus.phi() < phiMax_) &&
225  (Visible_Taus.pt() > ptMinMCTau_)) {
226  product_OneProng->push_back(Visible_Taus);
227  product_OneAndThreeProng->push_back(Visible_Taus);
228  product_Neutrina->push_back(Neutrino);
229  }
230  } else if (tauDecayMode == kThreeProng0pi0 || tauDecayMode == kThreeProng1pi0) {
231  if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
232  Visible_Taus.phi() < phiMax_) &&
233  (Visible_Taus.pt() > ptMinMCTau_)) {
234  product_ThreeProng->push_back(Visible_Taus);
235  product_OneAndThreeProng->push_back(Visible_Taus);
236  product_Neutrina->push_back(Neutrino);
237  }
238  } else if (tauDecayMode == kOther) {
239  if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
240  Visible_Taus.phi() < phiMax_) &&
241  (Visible_Taus.pt() > ptMinMCTau_)) {
242  product_Other->push_back(Visible_Taus);
243  }
244  }
245  }
246  }
247 
248  iEvent.put(std::move(product_Leptons), "LeptonicTauLeptons");
249  iEvent.put(std::move(product_Electrons), "LeptonicTauElectrons");
250  iEvent.put(std::move(product_Muons), "LeptonicTauMuons");
251  iEvent.put(std::move(product_OneProng), "HadronicTauOneProng");
252  iEvent.put(std::move(product_ThreeProng), "HadronicTauThreeProng");
253  iEvent.put(std::move(product_OneAndThreeProng), "HadronicTauOneAndThreeProng");
254  iEvent.put(std::move(product_Other), "TauOther");
255  iEvent.put(std::move(product_Neutrina), "Neutrina");
256  iEvent.put(std::move(product_MET), "MET");
257  iEvent.put(std::move(product_Mothers), "Mothers");
258 }
259 
260 // Helper Function
261 
264  int status,
265  int pdgId) const {
266  const GenParticleRefVector &daughterRefs = mother->daughterRefVector();
267 
268  for (GenParticleRefVector::const_iterator d = daughterRefs.begin(); d != daughterRefs.end(); ++d) {
269  if ((status == 0 || (*d)->status() == status) && (pdgId == 0 || std::abs((*d)->pdgId()) == pdgId)) {
270  products.push_back(*d);
271  } else
273  }
274 }
HiggsValidation_cfi.pdg_id
pdg_id
Definition: HiggsValidation_cfi.py:6
HLTTauMCProducer::getGenDecayProducts
void getGenDecayProducts(const reco::GenParticleRef &, reco::GenParticleRefVector &, int status=1, int pdgId=0) const
Definition: HLTTauMCProducer.cc:262
HLTTauMCProducer::ptMinMCMuon_
const double ptMinMCMuon_
Definition: HLTTauMCProducer.h:61
edm::StreamID
Definition: StreamID.h:30
HLTTauMCProducer::etaMax_
const double etaMax_
Definition: HLTTauMCProducer.h:63
HLTTauMCProducer::kOther
Definition: HLTTauMCProducer.h:53
genParticles2HepMC_cfi.genParticles
genParticles
Definition: genParticles2HepMC_cfi.py:4
HLTTauMCProducer::kOneProng0pi0
Definition: HLTTauMCProducer.h:48
HLTTauMCProducer::etaMin_
const double etaMin_
Definition: HLTTauMCProducer.h:63
HLTTauMCProducer::kElectron
Definition: HLTTauMCProducer.h:46
HLTTauMCProducer::MC_
const edm::EDGetTokenT< reco::GenParticleCollection > MC_
Definition: HLTTauMCProducer.h:57
CaloTowersParam_cfi.mc
mc
Definition: CaloTowersParam_cfi.py:8
mps_update.status
status
Definition: mps_update.py:69
edm
HLT enums.
Definition: AlignableModifier.h:19
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
edm::RefVector::begin
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223
HLTTauMCProducer::kThreeProng1pi0
Definition: HLTTauMCProducer.h:52
edm::RefVector< GenParticleCollection >
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
HLTTauMCProducer::phiMax_
const double phiMax_
Definition: HLTTauMCProducer.h:63
edm::Handle
Definition: AssociativeIterator.h:50
HLTTauMCProducer::produce
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
Definition: HLTTauMCProducer.cc:32
GenParticle
Definition: GenParticle.py:1
edm::es::products
ESProducts< std::remove_reference_t< TArgs >... > products(TArgs &&... args)
Definition: ESProducts.h:128
edm::Ref< GenParticleCollection >
reco::MET
Definition: MET.h:41
edm::RefVector::end
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
edm::RefVector::empty
bool empty() const
Is the RefVector empty.
Definition: RefVector.h:99
HLTTauMCProducer.h
HLTTauMCProducer::m_PDG_
const std::vector< int > m_PDG_
Definition: HLTTauMCProducer.h:62
HLTTauMCProducer::MCMET_
const edm::EDGetTokenT< reco::GenMETCollection > MCMET_
Definition: HLTTauMCProducer.h:58
edm::ParameterSet
Definition: ParameterSet.h:47
HLTTauMCProducer::ptMinMCTau_
const double ptMinMCTau_
Definition: HLTTauMCProducer.h:59
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:176
HLTTauMCProducer::kOneProng1pi0
Definition: HLTTauMCProducer.h:49
iEvent
int iEvent
Definition: GenABIO.cc:224
LorentzVectorCollection
std::vector< LorentzVector > LorentzVectorCollection
Definition: HLTTauMCProducer.h:31
EgammaValidation_cff.pdgId
pdgId
Definition: EgammaValidation_cff.py:118
edm::EventSetup
Definition: EventSetup.h:57
edm::RefVector::push_back
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
HepMCValidationHelper::genMet
TLorentzVector genMet(const HepMC::GenEvent *all, double etamin=-9999., double etamax=9999.)
Definition: HepMCValidationHelper.cc:278
reco::JetExtendedAssociation::LorentzVector
math::PtEtaPhiELorentzVectorF LorentzVector
Definition: JetExtendedAssociation.h:25
HLTTauMCProducer::phiMin_
const double phiMin_
Definition: HLTTauMCProducer.h:63
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
HLTTauMCProducer::kThreeProng0pi0
Definition: HLTTauMCProducer.h:51
LorentzVector
math::XYZTLorentzVector LorentzVector
Definition: HLTMuonMatchAndPlot.h:49
edm::RefVectorIterator
Definition: EDProductfwd.h:33
HLTTauMCProducer::ptMinMCElectron_
const double ptMinMCElectron_
Definition: HLTTauMCProducer.h:60
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
HLTTauMCProducer::HLTTauMCProducer
HLTTauMCProducer(const edm::ParameterSet &)
Definition: HLTTauMCProducer.cc:7
ztail.d
d
Definition: ztail.py:151
pi
const Double_t pi
Definition: trackSplitPlot.h:36
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::Event
Definition: Event.h:73
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
patZpeak.numMuons
numMuons
Definition: patZpeak.py:41
HLTTauDQMOffline_cfi.MET
MET
Definition: HLTTauDQMOffline_cfi.py:64
reco::GenParticle::isPromptDecayed
bool isPromptDecayed() const
Definition: GenParticle.h:55
edm::InputTag
Definition: InputTag.h:15
HLTTauMCProducer::kOneProng2pi0
Definition: HLTTauMCProducer.h:50
HLTTauMCProducer::kMuon
Definition: HLTTauMCProducer.h:47