CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
FlavorHistoryProducer Class Reference

#include <FlavorHistoryProducer.h>

Inheritance diagram for FlavorHistoryProducer:
edm::EDProducer edm::ProducerBase edm::ProductRegistryHelper

Public Member Functions

 FlavorHistoryProducer (const edm::ParameterSet &)
 constructor More...
 
 ~FlavorHistoryProducer ()
 destructor More...
 
- Public Member Functions inherited from edm::EDProducer
 EDProducer ()
 
virtual ~EDProducer ()
 
- Public Member Functions inherited from edm::ProducerBase
 ProducerBase ()
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
boost::function< void(const
BranchDescription &)> 
registrationCallback () const
 used by the fwk to register list of products More...
 
virtual ~ProducerBase ()
 

Private Member Functions

void getAncestors (const reco::Candidate &c, std::vector< reco::Candidate const * > &moms)
 
reco::CandidateView::const_iterator getClosestJet (edm::Handle< reco::CandidateView > const &pJets, reco::Candidate const &parton) const
 
reco::Candidate const * getSister (const reco::Candidate &c)
 
void produce (edm::Event &e, const edm::EventSetup &)
 process one event More...
 

Private Attributes

double etaMaxParticle_
 
double etaMaxShower_
 
std::string flavorHistoryName_
 
double matchDR_
 
edm::InputTag matchedSrc_
 
int pdgIdToSelect_
 
double ptMinParticle_
 
double ptMinShower_
 
edm::InputTag src_
 
bool verbose_
 

Additional Inherited Members

- Public Types inherited from edm::EDProducer
typedef EDProducer ModuleType
 
typedef WorkerT< EDProducerWorkerType
 
- Public Types inherited from edm::ProducerBase
typedef
ProductRegistryHelper::TypeLabelList 
TypeLabelList
 
- Static Public Member Functions inherited from edm::EDProducer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &descriptions)
 
- Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const * currentContext () const
 
- Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)
 

Detailed Description

class

Author
Stephen Mrenna, FNAL
Version
$Id: FlavorHistoryProducer.cc,v 1.0

Definition at line 57 of file FlavorHistoryProducer.h.

Constructor & Destructor Documentation

FlavorHistoryProducer::FlavorHistoryProducer ( const edm::ParameterSet p)

constructor

Definition at line 38 of file FlavorHistoryProducer.cc.

References flavorHistoryName_.

38  :
39  src_( p.getParameter<InputTag>( "src" ) ),
40  matchedSrc_( p.getParameter<InputTag>( "matchedSrc") ),
41  matchDR_ ( p.getParameter<double> ("matchDR") ),
42  pdgIdToSelect_( p.getParameter<int> ("pdgIdToSelect") ),
43  ptMinParticle_( p.getParameter<double>( "ptMinParticle") ),
44  ptMinShower_( p.getParameter<double>( "ptMinShower") ),
45  etaMaxParticle_( p.getParameter<double>( "etaMaxParticle" )),
46  etaMaxShower_( p.getParameter<double>( "etaMaxShower" )),
47  flavorHistoryName_( p.getParameter<string>("flavorHistoryName") ),
48  verbose_( p.getUntrackedParameter<bool>( "verbose" ) )
49 {
50  produces<FlavorHistoryEvent >(flavorHistoryName_);
51 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
FlavorHistoryProducer::~FlavorHistoryProducer ( )

destructor

Definition at line 53 of file FlavorHistoryProducer.cc.

53  {
54 }

Member Function Documentation

void FlavorHistoryProducer::getAncestors ( const reco::Candidate c,
std::vector< reco::Candidate const * > &  moms 
)
private

Definition at line 329 of file FlavorHistoryProducer.cc.

References trackerHits::c, reco::Candidate::mother(), and reco::Candidate::numberOfMothers().

Referenced by produce().

331 {
332 
333  if( c.numberOfMothers() == 1 ) {
334  const Candidate * dau = &c;
335  const Candidate * mom = c.mother();
336  while ( dau->numberOfMothers() != 0) {
337  moms.push_back( dau );
338  dau = mom ;
339  mom = dau->mother();
340  }
341  }
342 }
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
CandidateView::const_iterator FlavorHistoryProducer::getClosestJet ( edm::Handle< reco::CandidateView > const &  pJets,
reco::Candidate const &  parton 
) const
private

Definition at line 346 of file FlavorHistoryProducer.cc.

References deltaR(), j, matchDR_, and reco::Candidate::p4().

Referenced by produce().

348 {
349  double dr = matchDR_;
350  CandidateView::const_iterator j = pJets->begin(),
351  jend = pJets->end();
352  CandidateView::const_iterator bestJet = pJets->end();
353  for ( ; j != jend; ++j ) {
354  double dri = deltaR( parton.p4(), j->p4() );
355  if ( dri < dr ) {
356  dr = dri;
357  bestJet = j;
358  }
359  }
360  return bestJet;
361 }
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
int j
Definition: DBlmapReader.cc:9
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
reco::Candidate const* FlavorHistoryProducer::getSister ( const reco::Candidate c)
private
void FlavorHistoryProducer::produce ( edm::Event e,
const edm::EventSetup  
)
privatevirtual

process one event

Implements edm::EDProducer.

Definition at line 56 of file FlavorHistoryProducer.cc.

References abs, gather_cfg::cout, eta(), etaMaxShower_, spr::find(), flavorHistoryName_, newFWLiteAna::found, getAncestors(), edm::Event::getByLabel(), getClosestJet(), i, j, findQualityFiles::jj, matchedSrc_, reco::Candidate::numberOfMothers(), AlCaHLTBitMon_ParallelJobs::p, benchmark_cfg::pdgId, reco::Candidate::pdgId(), pdgIdToSelect_, ptMinShower_, edm::Event::put(), src_, reco::Candidate::status(), ntuplemaker::status, and verbose_.

57 {
58  if ( verbose_ ) cout << "---------- Hello from FlavorHistoryProducer! -----" << endl;
59 
60  // Get a handle to the particle collection (OwnVector)
61  Handle<CandidateView > particlesViewH;
62  evt.getByLabel( src_, particlesViewH );
63 
64  Handle<CandidateView> matchedView;
65  evt.getByLabel( matchedSrc_, matchedView );
66 
67  // Copy the View to an vector for easier iterator manipulation convenience
68  vector<const Candidate* > particles;
69  for( CandidateView::const_iterator p = particlesViewH->begin(); p != particlesViewH->end(); ++p ) {
70  particles.push_back(&*p);
71  }
72 
73  // List of indices for the partons to add
74  vector<int> partonIndices;
75  // List of progenitors for those partons
76  vector<int> progenitorIndices;
77  // List of sisters for those partons
78  vector<int> sisterIndices;
79  // Flavor sources
80  vector<FlavorHistory::FLAVOR_T> flavorSources;
81 
82  // Make a new flavor history vector
83  auto_ptr<FlavorHistoryEvent > flavorHistoryEvent ( new FlavorHistoryEvent () ) ;
84 
85  // ------------------------------------------------------------
86  // Loop over partons
87  // ------------------------------------------------------------
89  vector<const Candidate* >::size_type j_max = particles.size();
90  for( j=0; j<j_max; ++j ) {
91 
92  // Get the candidate
93  const Candidate *p = particles[j];
94  // Set up indices that we'll need for the flavor history
96  vector<Candidate const *>::size_type progenitorIndex=j_max;
97  bool foundProgenitor = false;
98  FlavorHistory::FLAVOR_T flavorSource=FlavorHistory::FLAVOR_NULL;
99 
100 
101  int idabs = std::abs( (p)->pdgId() );
102  int nDa = (p)->numberOfDaughters();
103 
104  // Check if we have a status 2 or 3 particle, which is a parton before the string.
105  // Only consider quarks.
106  if ( idabs <= FlavorHistory::tQuarkId && p->status() == 2 &&
107  idabs == pdgIdToSelect_ ) {
108  // Ensure the parton in question has daughters
109  if ( nDa > 0 ) {
110  // Ensure the parton passes some minimum kinematic cuts
111  if((p)->pt() > ptMinShower_ && fabs((p)->eta())<etaMaxShower_) {
112 
113  if(verbose_) cout << "--------------------------" << endl;
114  if(verbose_) cout << "Processing particle " << j << ", status = " << p->status() << ", pdg id = " << p->pdgId() << endl;
115 
116 
117  //Main (but simple) workhorse to get all ancestors
118  vector<Candidate const *> allParents;
119  getAncestors( *p, allParents );
120 
121  if(verbose_) {
122  cout << "Parents : " << endl;
123  vector<Candidate const *>::const_iterator iprint = allParents.begin(),
124  iprintend = allParents.end();
125  for ( ; iprint != iprintend; ++iprint )
126  cout << " status = " << (*iprint)->status() << ", pdg id = " << (*iprint)->pdgId() << ", pt = " << (*iprint)->pt() << endl;
127  }
128 
129  // -------------------------------------------------------------
130  // Identify the ancestry of the Quark
131  //
132  //
133  // Matrix Element:
134  // Status 3 parent with precisely 2 "grandparents" that
135  // is outside of the "initial" section (0-5) that has the
136  // same ID as the status 2 parton in question.
137  //
138  // Flavor excitation:
139  // If we find only one outgoing parton.
140  //
141  // Gluon splitting:
142  // Parent is a quark of a different flavor than the parton
143  // in question, or a gluon.
144  // Can come from either ISR or FSR.
145  //
146  // True decay:
147  // Decays from a resonance like top, Higgs, etc.
148  // -------------------------------------------------------------
149  vector<Candidate const *>::size_type a_size = allParents.size();
150 
151  //
152  // Loop over all the ancestors of this parton and find the progenitor.
153  //
154  for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
155  const Candidate * aParent=allParents[i];
156  vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
157 
158 
159  // Get the index of the progenitor candidate
160  progenitorIndex = found - particles.begin();
161 
162  int aParentId = std::abs(aParent->pdgId());
163 
164  // -----------------------------------------------------------------------
165  // Here we examine particles that were produced after the collision
166  // -----------------------------------------------------------------------
167  if( progenitorIndex > 5 ) {
168  // Here is where we have a matrix element process.
169  // The signature is to have a status 3 parent with precisely 2 parents
170  // that is outside of the "initial" section (0-5) that has the
171  // same ID as the status 2 parton in question.
172  // NOTE: If this parton has no sister, then this will be classified as
173  // a flavor excitation. Often the only difference is that the matrix element
174  // cases have a sister, while the flavor excitation cases do not.
175  // If we do not find a sister below, this will be classified as a flavor
176  // excitation.
177  if( aParent->numberOfMothers() == 2 &&
178  aParent->pdgId() == p->pdgId() && aParent->status() == 3 ) {
179  if(verbose_) cout << "Matrix Element progenitor" << endl;
180  if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_ME;
181  foundProgenitor = true;
182  }
183  // Here we have a true decay. Parent is not a quark or a gluon.
184  else if( (aParentId>pdgIdToSelect_ && aParentId<FlavorHistory::gluonId) ||
185  aParentId > FlavorHistory::gluonId ) {
186  if(verbose_) cout << "Flavor decay progenitor" << endl;
187  if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_DECAY;
188  foundProgenitor = true;
189  }
190  }// end if progenitorIndex > 5
191  }// end loop over parents
192 
193  // Otherwise, classify it as gluon splitting. Take the first status 3 parton with 1 parent
194  // that you see as the progenitor
195  if ( !foundProgenitor ) {
196  if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_GS;
197  // Now get the parton with only one parent (the proton) and that is the progenitor
198  for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
199  const Candidate * aParent=allParents[i];
200  vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
201  // Get the index of the progenitor candidate
202  progenitorIndex = found - particles.begin();
203 
204  if ( aParent->numberOfMothers() == 1 && aParent->status() == 3 ) {
205  foundProgenitor = true;
206  }
207  }// end loop over parents
208  }// end if we haven't found a progenitor, and if we haven't found a status 3 parton of the same flavor
209  // (i.e. end if examining gluon splitting)
210  }// End if this parton passes some minimal kinematic cuts
211  }// End if this particle is status 2 (has strings as daughters)
212 
213 
214 
215  // Make sure we've actually found a progenitor
216  if ( !foundProgenitor ) progenitorIndex = j_max;
217 
218  // We've found the particle, add to the list
219 
220  partonIndices.push_back( partonIndex );
221  progenitorIndices.push_back( progenitorIndex );
222  flavorSources.push_back(flavorSource);
223  sisterIndices.push_back( -1 ); // set below
224 
225  }// End if this particle was a status==2 parton
226  }// end loop over particles
227 
228  // Now add sisters.
229  // Also if the event is preliminarily classified as "matrix element", check to
230  // make sure that they have a sister. If not, it is flavor excitation.
231 
232  if ( verbose_ ) cout << "Making sisters" << endl;
233  // First make sure nothing went terribly wrong:
234  if ( partonIndices.size() == progenitorIndices.size() && partonIndices.size() > 0 ) {
235  // Now loop over the candidates
236  for ( unsigned int ii = 0; ii < partonIndices.size(); ++ii ) {
237  // Get the iith particle
238  const Candidate * candi = particles[partonIndices[ii]];
239  if ( verbose_ ) cout << " Sister candidate particle 1: " << ii << ", pdgid = " << candi->pdgId() << ", status = " << candi->status() << endl;
240  // Get the iith progenitor
241  // Now loop over the other flavor history candidates and
242  // attempt to find a sister to this one
243  for ( unsigned int jj = 0; jj < partonIndices.size(); ++jj ) {
244  if ( ii != jj ) {
245  const Candidate * candj = particles[partonIndices[jj]];
246  if ( verbose_ ) cout << " Sister candidate particle 2: " << jj << ", pdgid = " << candj->pdgId() << ", status = " << candj->status() << endl;
247  // They should be opposite in pdgid and have the same status, and
248  // the same progenitory.
249  if ( candi->pdgId() == -1 * candj->pdgId() && candi->status() == candj->status() ) {
250  sisterIndices[ii] = partonIndices[jj];
251  if ( verbose_ ) cout << "Parton " << partonIndices[ii] << " has sister " << sisterIndices[ii] << endl;
252  }
253  }
254  }
255 
256  // Here, ensure that there is a sister. Otherwise this is "flavor excitation"
257  if ( sisterIndices[ii] < 0 ) {
258  if ( verbose_ ) cout << "No sister. Classified as flavor excitation" << endl;
259  flavorSources[ii] = FlavorHistory::FLAVOR_EXC;
260  }
261 
262  if ( verbose_ ) cout << "Getting matched jet" << endl;
263  // Get the closest match in the matched object collection
264  CandidateView::const_iterator matched = getClosestJet( matchedView, *candi );
265  CandidateView::const_iterator matchedBegin = matchedView->begin();
266  CandidateView::const_iterator matchedEnd = matchedView->end();
267 
268 
269  if ( verbose_ ) cout << "Getting sister jet" << endl;
270  // Get the closest sister in the matched object collection, if sister is found
272  ( sisterIndices[ii] >= 0 && static_cast<unsigned int>(sisterIndices[ii]) < particles.size() ) ?
273  getClosestJet( matchedView, *particles[sisterIndices[ii]] ) :
274  matchedEnd ;
275 
276  if ( verbose_ ) cout << "Making matched shallow clones : ";
277  ShallowClonePtrCandidate matchedCand ;
278  if ( matched != matchedEnd ) {
279  if ( verbose_ ) cout << " found" << endl;
280  matchedCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, matched - matchedBegin ) );
281  } else {
282  if ( verbose_ ) cout << " NOT found" << endl;
283  }
284 
285  if ( verbose_ ) cout << "Making sister shallow clones : ";
286  ShallowClonePtrCandidate sisterCand;
287  if ( sister != matchedEnd ) {
288  if ( verbose_ ) cout << " found" << endl;
289  sisterCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, sister - matchedBegin ) );
290  } else {
291  if ( verbose_ ) cout << " NOT found" << endl;
292  }
293 
294  if ( verbose_ ) cout << "Making history object" << endl;
295  // Now create the flavor history object
296  FlavorHistory history (flavorSources[ii],
297  particlesViewH,
298  partonIndices[ii], progenitorIndices[ii], sisterIndices[ii],
299  matchedCand,
300  sisterCand );
301  if ( verbose_ ) cout << "Adding flavor history : " << history << endl;
302  flavorHistoryEvent->push_back( history );
303  }
304  }
305 
306  // If we got any partons, cache them and then write them out
307  if ( flavorHistoryEvent->size() > 0 ) {
308 
309  // Calculate some nice variables for FlavorHistoryEvent
310  if ( verbose_ ) cout << "Caching flavor history event" << endl;
311  flavorHistoryEvent->cache();
312 
313  if ( verbose_ ) {
314  cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryEvent->size() << endl;
315  vector<FlavorHistory>::const_iterator i = flavorHistoryEvent->begin(),
316  iend = flavorHistoryEvent->end();
317  for ( ; i !=iend; ++i ) {
318  cout << *i << endl;
319  }
320  }
321  }
322 
323  // Now add the flavor history to the event record
324  evt.put( flavorHistoryEvent, flavorHistoryName_ );
325 }
int i
Definition: DBlmapReader.cc:9
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
virtual int status() const =0
status word
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
#define abs(x)
Definition: mlp_lapack.h:159
T eta() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
uint16_t size_type
reco::CandidateView::const_iterator getClosestJet(edm::Handle< reco::CandidateView > const &pJets, reco::Candidate const &parton) const
int j
Definition: DBlmapReader.cc:9
virtual int pdgId() const =0
PDG identifier.
edm::Ptr< Candidate > CandidatePtr
persistent reference to an object in a collection of Candidate objects
Definition: CandidateFwd.h:25
void getAncestors(const reco::Candidate &c, std::vector< reco::Candidate const * > &moms)
tuple cout
Definition: gather_cfg.py:121
tuple status
Definition: ntuplemaker.py:245

Member Data Documentation

double FlavorHistoryProducer::etaMaxParticle_
private

Definition at line 84 of file FlavorHistoryProducer.h.

double FlavorHistoryProducer::etaMaxShower_
private

Definition at line 85 of file FlavorHistoryProducer.h.

Referenced by produce().

std::string FlavorHistoryProducer::flavorHistoryName_
private

Definition at line 86 of file FlavorHistoryProducer.h.

Referenced by FlavorHistoryProducer(), and produce().

double FlavorHistoryProducer::matchDR_
private

Definition at line 80 of file FlavorHistoryProducer.h.

Referenced by getClosestJet().

edm::InputTag FlavorHistoryProducer::matchedSrc_
private

Definition at line 79 of file FlavorHistoryProducer.h.

Referenced by produce().

int FlavorHistoryProducer::pdgIdToSelect_
private

Definition at line 81 of file FlavorHistoryProducer.h.

Referenced by produce().

double FlavorHistoryProducer::ptMinParticle_
private

Definition at line 82 of file FlavorHistoryProducer.h.

double FlavorHistoryProducer::ptMinShower_
private

Definition at line 83 of file FlavorHistoryProducer.h.

Referenced by produce().

edm::InputTag FlavorHistoryProducer::src_
private

Definition at line 78 of file FlavorHistoryProducer.h.

Referenced by produce().

bool FlavorHistoryProducer::verbose_
private

Definition at line 87 of file FlavorHistoryProducer.h.

Referenced by produce().