---

FlavorHistoryProducer Class Reference

#include <PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h>

Inheritance diagram for FlavorHistoryProducer:

List of all members.

Public Member Functions
	FlavorHistoryProducer (const edm::ParameterSet &)
	constructor
	~FlavorHistoryProducer ()
	destructor
Private Member Functions
void	beginJob (const edm::EventSetup &)
	module init at begin of job
void	getAncestors (const reco::Candidate &c, std::vector< reco::Candidate const * > &moms)
void	produce (edm::Event &e, const edm::EventSetup &)
	process one event
Private Attributes
double	etaMaxParticle_
double	etaMaxShower_
std::string	flavorHistoryName_
int	pdgIdToSelect_
double	ptMinParticle_
double	ptMinShower_
edm::InputTag	src_
bool	verbose_

---

Detailed Description

Definition at line 64 of file FlavorHistoryProducer.h.

---

Constructor & Destructor Documentation

FlavorHistoryProducer::FlavorHistoryProducer

(

const edm::ParameterSet &

p

)

constructor

Definition at line 74 of file FlavorHistoryProducer.cc.

References flavorHistoryName_.

                                                                     :
  src_( p.getParameter<InputTag>( "src" ) ),
  pdgIdToSelect_( p.getParameter<int> ("pdgIdToSelect") ),
  ptMinParticle_( p.getParameter<double>( "ptMinParticle") ),  
  ptMinShower_( p.getParameter<double>( "ptMinShower") ),  
  etaMaxParticle_( p.getParameter<double>( "etaMaxParticle" )),  
  etaMaxShower_( p.getParameter<double>( "etaMaxShower" )),
  flavorHistoryName_( p.getParameter<string>("flavorHistoryName") ),
  verbose_( p.getUntrackedParameter<bool>( "verbose" ) )
{
  produces<vector<FlavorHistory> >(flavorHistoryName_);
}

FlavorHistoryProducer::~FlavorHistoryProducer

(

)

destructor

Definition at line 87 of file FlavorHistoryProducer.cc.

                                              { 
}

---

Member Function Documentation

void FlavorHistoryProducer::beginJob

(

const edm::EventSetup &

es

)

[private, virtual]

module init at begin of job

Reimplemented from edm::EDProducer.

Definition at line 90 of file FlavorHistoryProducer.cc.

                                                            {
}

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

Referenced by produce().

void FlavorHistoryProducer::produce	(	edm::Event &	e,
		const edm::EventSetup &
	)			[private, virtual]

process one event

Implements edm::EDProducer.

Definition at line 93 of file FlavorHistoryProducer.cc.

References funct::abs(), GenMuonPlsPt100GeV_cfg::cout, reco::Candidate::daughter(), lat::endl(), eta, etaMaxShower_, find(), flavorHistoryName_, getAncestors(), edm::Event::getByLabel(), i, j, reco::Candidate::mother(), reco::Candidate::numberOfDaughters(), reco::Candidate::numberOfMothers(), p, reco::Particle::pdgId(), pdgIdToSelect_, ptMinShower_, edm::Event::put(), src_, StDecayID::status, reco::Particle::status(), and verbose_.

{

  if ( verbose_ ) cout << "Producing flavor history" << endl;
  
  // Get a handle to the particle collection (OwnVector)
  Handle<CandidateView > particlesViewH;
  evt.getByLabel( src_, particlesViewH );

  // Copy the View to an vector for easier iterator manipulation convenience
  vector<const Candidate* > particles;
  for( CandidateView::const_iterator p = particlesViewH->begin();  p != particlesViewH->end(); ++p ) {
    particles.push_back(&*p);
  }
  
  // Make a new flavor history vector
  auto_ptr<vector<FlavorHistory> > flavorHistoryVector ( new vector<FlavorHistory> () ) ;

  // ------------------------------------------------------------
  // Loop over partons
  // ------------------------------------------------------------
  vector<const Candidate* >::size_type j;
  vector<const Candidate* >::size_type j_max = particles.size();
  for( j=0; j<j_max; ++j ) {

    if ( verbose_ ) cout << "Processing particle " << j << endl;

    // Get the candidate
    const Candidate *p = particles[j];
    // Set up indices that we'll need for the flavor history
    vector<Candidate const *>::size_type partonIndex=j;
    vector<Candidate const *>::size_type progenitorIndex=0;
    vector<Candidate const *>::size_type sisterIndex=0;
    bool foundProgenitor = false; 
    bool foundSister = false;
    FlavorHistory::FLAVOR_T flavorSource=FlavorHistory::FLAVOR_NULL;


    int idabs = abs( (p)->pdgId() );
    int nDa = (p)->numberOfDaughters();

    // Check if we have a status 2 or 3 particle, which is a parton before the string.
    // Only consider quarks. 
    if ( idabs <= FlavorHistory::tQuarkId && p->status() == 2 ) {
      // Ensure the parton in question has daughters 
      if ( nDa > 0 ) {
        // Ensure the parton passes some minimum kinematic cuts
        if((p)->pt() > ptMinShower_ && fabs((p)->eta())<etaMaxShower_) {


          if(verbose_) cout << "--------------------------" << endl;
          if(verbose_) cout << "Processing particle " << j  << " = " << *p << endl;


          //Main (but simple) workhorse to get all ancestors
          vector<Candidate const *> allParents;
          getAncestors( *p, allParents );
            
          if(verbose_) {
            cout << "Parents : " << endl;
            vector<Candidate const *>::const_iterator iprint = allParents.begin(),
              iprintend = allParents.end();
            for ( ; iprint != iprintend; ++iprint ) 
              cout << **iprint << endl;
          }
          

          // -------------------------------------------------------------
          // Identify the ancestry of the Quark
          // 
          // 
          // Matrix Element:
          //    Status 3 parent with precisely 2 "grandparents" that
          //    is outside of the "initial" section (0-5) that has the
          //    same ID as the status 2 parton in question. 
          //    NOTE: This is not the actual ultimate progenitor,
          //    but this is the signature of matrix element decays.
          //    The ultimate progenitor is the parent of the status 3
          //    parton.
          //
          // Flavor excitation:
          //    Almost the same as the matrix element classification,
          //    but has only one outgoing parton product instead of two.
          //
          // Gluon splitting:
          //    Parent is a quark of a different flavor than the parton
          //    in question, or a gluon. Can come from either ISR or FSR.
          //
          // True decay:
          //    Decays from a resonance like top, Higgs, etc.
          // -------------------------------------------------------------
          vector<Candidate const *>::size_type a_size = allParents.size();
          int parentIndex=0;

          // 
          // Loop over all the ancestors of this parton and find the progenitor.
          // 
          for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i,++parentIndex ) {
            const Candidate * aParent=allParents[i];
            vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);


            // Get the index of the progenitor candidate
            progenitorIndex = found - particles.begin();

            int aParentId = abs(aParent->pdgId());

            // -----------------------------------------------------------------------
            // Here we examine particles that were produced after the collision
            // -----------------------------------------------------------------------
            if( aParent->numberOfMothers() == 2 && progenitorIndex > 5 ) {
              // Here is where we have a matrix element process.
              // The signature is to have a status 3 parent with precisely 2 parents
              // that is outside of the "initial" section (0-5) that has the
              // same ID as the status 2 parton in question.
              // NOTE: This is not the actual ultimate progenitor, but this is
              // the signature of matrix element decays. The ultimate progentitor
              // is the parent of the status 3 parton. 
              // ALSO NOTE: If this parton has no sister, then this will be classified as
              // a flavor excitation. The only difference, since the initial states are
              // mostly gluons, is that the matrix element cases have a sister,
              // while the flavor excitation cases do not. 
              // If we do not find a sister below, this will be classified as a flavor
              // excitation. 
              if( aParentId == pdgIdToSelect_ ) {
                if(verbose_) cout << "Matrix Element progenitor" << endl;
                flavorSource = FlavorHistory::FLAVOR_ME;

                // The "true" progenitor is the next parent in the list (the parent of this
                // progenitor).
                if ( i != a_size - 1 ) {
                  const Candidate * progenitorCand = allParents[i+1];
                  vector<Candidate const *>::const_iterator foundAgain = find( particles.begin(),
                                                                               particles.end(),
                                                                               progenitorCand );
                  progenitorIndex = foundAgain - particles.begin();
                  foundProgenitor = true;

                } else {
                  edm::LogWarning("UnexpectedFormat") << "Error! Parentage information in FlavorHistoryProducer is not what is expected";
                  
                  cout << "Particle : " << *p << endl;
                  cout << "Parents : " << endl;
                  vector<Candidate const *>::const_iterator iprint = allParents.begin(),
                    iprintend = allParents.end();
                  for ( ; iprint != iprintend; ++iprint ) 
                    cout << **iprint << endl;

                  foundProgenitor = false;

                }
              } 
              // Here we have a gluon splitting from final state radiation. 
              // The parent is a quark of a different flavor, or a gluon, in the
              // final state. 
              else if( (aParentId > 0 && aParentId < FlavorHistory::tQuarkId ) || aParentId==FlavorHistory::gluonId ) {
                if(verbose_) cout << "Gluon splitting progenitor" << endl;
                flavorSource = FlavorHistory::FLAVOR_GS;
                foundProgenitor = true;
              }
              // Here we have a true decay. Parent is not a quark or a gluon.
              else if( (aParentId>pdgIdToSelect_ && aParentId<FlavorHistory::gluonId) || aParentId > FlavorHistory::gluonId ) {
                if(verbose_) cout << "Flavor decay progenitor" << endl;
                flavorSource = FlavorHistory::FLAVOR_DECAY;
                foundProgenitor = true;
              }
            }

            // -----------------------------------------------------------------------
            // Here we examine particles that were produced before the collision
            // -----------------------------------------------------------------------
            else if( progenitorIndex <= 5 ) {
              // Parent has a quark daughter equal and opposite to this: ISR
              if( aParent->numberOfDaughters() > 0 && 
                  aParent->daughter(0)->pdgId() == -1 * p->pdgId()  ) {
                if(verbose_) cout << "Gluon splitting progenitor" << endl;
                flavorSource = FlavorHistory::FLAVOR_GS;
              }
              // Otherwise, this is flavor excitation. Rarely happens because
              // mostly the initial state is gluons which will be caught by the
              // "matrix element" version above. 
              else {              
                if(verbose_) cout << "Flavor excitation progenitor" << endl;
                flavorSource = FlavorHistory::FLAVOR_EXC;
              }
              foundProgenitor = true;
            }
          }// End loop over all parents of this parton to find progenitor

            

          // 
          // Now find sister of this particle if there is one
          // 
          if ( foundProgenitor && progenitorIndex >= 0 ) {
            // Get the progenitor
            const Candidate * progenitorCand = particles[progenitorIndex];

            if ( verbose_ ) cout << "Found progenitor: " << *progenitorCand << endl;

            // Here is the normal case of a sister around
            if ( progenitorCand->numberOfDaughters() >= 2 ) {
              const Candidate * sisterCand = 0;
              
              for ( unsigned int iida = 0; iida < progenitorCand->numberOfDaughters(); ++iida ) {
                const Candidate * dai = progenitorCand->daughter(iida);

                if ( verbose_ ) cout << "Sister candidate " << *dai << endl;
                
                if ( dai->pdgId() == -1 * p->pdgId() ) {
                  if ( verbose_ ) cout << "Found actual sister" << endl;
                  sisterCand = dai;
                  foundSister = true;
                }
              }
                
              if ( foundSister ) {
                // Find index of daughter in master list
                vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),sisterCand);
                sisterIndex = found - particles.begin();
                if(verbose_) cout << "Sister index = " << sisterIndex << endl;
                if ( found != particles.end() )
                  if(verbose_) cout << "Sister = " << **found << endl;
              } // end if found sister
            }
            // Here is if we have a "transient" decay in the code that isn't
            // really a decay, so we need to look at the parent of the progenitor
            else {
              const Candidate * grandProgenitorCand = progenitorCand->mother(0);
              const Candidate * sisterCand = 0;

              if ( verbose_ ) cout << "Looking for sister, progenitor is " << *progenitorCand << endl;
            
              // Make sure the progenitor has two daughters
              if ( grandProgenitorCand->numberOfDaughters() >= 2 ) {

                for ( unsigned int iida = 0; iida < grandProgenitorCand->numberOfDaughters(); ++iida ) {
                  const Candidate * dai = grandProgenitorCand->daughter(iida);

                  if ( verbose_ ) cout << "Looking for sister " << *dai << endl;
                
                  if ( dai->pdgId() == -1 * p->pdgId() ) {
                    if ( verbose_ ) cout << "Found sister" << endl;
                    sisterCand = dai;
                    foundSister = true;
                  }
                }
                
                if ( foundSister ) {
                  // Find index of daughter in master list
                  vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),sisterCand);
                  sisterIndex = found - particles.begin();
                  if(verbose_) cout << "Sister index = " << sisterIndex << endl;
                  if ( found != particles.end() )
                    if(verbose_) cout << "Sister = " << **found << endl;
                } // end if found sister
              } // End of have at least 2 grand progenitor daughters
            } // End if we have to look at parents of progenitor to find sister

          } // end if found progenitor

          // ------
          // Here, we change the type from matrix element to flavor excitation
          // if there are no sisters present. 
          // ------
          if ( flavorSource == FlavorHistory::FLAVOR_ME && !foundSister ) {
            flavorSource = FlavorHistory::FLAVOR_EXC;
          }
          
        }// End if this parton passes some minimal kinematic cuts
      }// End if this particle has strings as daughters
    }// End if this particle was a status==2 parton

    // Make sure we've actually found a sister and a progenitor
    if ( !foundProgenitor ) progenitorIndex = 0;
    if ( !foundSister ) sisterIndex = 0;

    // We've found the particle, add to the list (status 2 only)
    if ( idabs == pdgIdToSelect_ && p->status() == 2 ) 
      flavorHistoryVector->push_back( FlavorHistory( flavorSource, particlesViewH, partonIndex, progenitorIndex, sisterIndex ) ); 
  }


//   ValueMap<FlavorHistory>::Filler filler(*flavorHistory);
//   filler.insert( particlesViewH, flavorHistoryVector.begin(), flavorHistoryVector.end()  );
//   filler.fill();
  // Now add the flavor history to the event record
  if ( verbose_ ) {
    cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryVector->size() << endl;
    vector<FlavorHistory>::const_iterator i = flavorHistoryVector->begin(),
      iend = flavorHistoryVector->end();
    for ( ; i !=iend; ++i ) {
      cout << *i << endl;
    }
  }
  evt.put( flavorHistoryVector, flavorHistoryName_ );
}

---

Member Data Documentation

double FlavorHistoryProducer::etaMaxParticle_ [private]

Definition at line 86 of file FlavorHistoryProducer.h.

double FlavorHistoryProducer::etaMaxShower_ [private]

Definition at line 87 of file FlavorHistoryProducer.h.

Referenced by produce().

std::string FlavorHistoryProducer::flavorHistoryName_ [private]

Definition at line 88 of file FlavorHistoryProducer.h.

Referenced by FlavorHistoryProducer(), and produce().

int FlavorHistoryProducer::pdgIdToSelect_ [private]

Definition at line 83 of file FlavorHistoryProducer.h.

Referenced by produce().

double FlavorHistoryProducer::ptMinParticle_ [private]

Definition at line 84 of file FlavorHistoryProducer.h.

double FlavorHistoryProducer::ptMinShower_ [private]

Definition at line 85 of file FlavorHistoryProducer.h.

Referenced by produce().

edm::InputTag FlavorHistoryProducer::src_ [private]

Definition at line 82 of file FlavorHistoryProducer.h.

Referenced by produce().

bool FlavorHistoryProducer::verbose_ [private]

Definition at line 89 of file FlavorHistoryProducer.h.

Referenced by produce().

---

The documentation for this class was generated from the following files:

• PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h
• PhysicsTools/HepMCCandAlgos/plugins/FlavorHistoryProducer.cc

---