CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/PhysicsTools/HepMCCandAlgos/plugins/FlavorHistoryProducer.cc

Go to the documentation of this file.
00001 // This file was removed but it should not have been.
00002 // This comment is to restore it. 
00003 
00004 #include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h"
00005 
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "DataFormats/Math/interface/deltaR.h"
00008 
00009 // #include "DataFormats/Common/interface/ValueMap.h"
00010 // #include <iterators>
00011 
00012 
00013 // -------------------------------------------------------------
00014 // Identify the ancestry of the Quark
00015 // 
00016 // 
00017 // Matrix Element:
00018 //    Status 3 parent with precisely 2 "grandparents" that
00019 //    is outside of the "initial" section (0-5) that has the
00020 //    same ID as the status 2 parton in question. 
00021 //
00022 // Flavor excitation:
00023 //    If we find only one outgoing parton.
00024 //
00025 // Gluon splitting:
00026 //    Parent is a quark of a different flavor than the parton
00027 //    in question, or a gluon. 
00028 //    Can come from either ISR or FSR.
00029 //
00030 // True decay:
00031 //    Decays from a resonance like top, Higgs, etc.
00032 // -------------------------------------------------------------
00033 
00034 using namespace std;
00035 using namespace reco;
00036 using namespace edm;
00037 
00038 FlavorHistoryProducer::FlavorHistoryProducer( const ParameterSet & p ) :
00039   src_( p.getParameter<InputTag>( "src" ) ),
00040   matchedSrc_( p.getParameter<InputTag>( "matchedSrc") ),
00041   matchDR_ ( p.getParameter<double> ("matchDR") ),
00042   pdgIdToSelect_( p.getParameter<int> ("pdgIdToSelect") ),
00043   ptMinParticle_( p.getParameter<double>( "ptMinParticle") ),  
00044   ptMinShower_( p.getParameter<double>( "ptMinShower") ),  
00045   etaMaxParticle_( p.getParameter<double>( "etaMaxParticle" )),  
00046   etaMaxShower_( p.getParameter<double>( "etaMaxShower" )),
00047   flavorHistoryName_( p.getParameter<string>("flavorHistoryName") ),
00048   verbose_( p.getUntrackedParameter<bool>( "verbose" ) )
00049 {
00050   produces<FlavorHistoryEvent >(flavorHistoryName_);
00051 }
00052 
00053 FlavorHistoryProducer::~FlavorHistoryProducer() { 
00054 }
00055 
00056 void FlavorHistoryProducer::produce( Event& evt, const EventSetup& ) 
00057 {
00058   if ( verbose_ ) cout << "---------- Hello from FlavorHistoryProducer! -----" << endl;
00059   
00060   // Get a handle to the particle collection (OwnVector)
00061   Handle<CandidateView > particlesViewH;
00062   evt.getByLabel( src_, particlesViewH );
00063 
00064   Handle<CandidateView> matchedView;
00065   evt.getByLabel( matchedSrc_, matchedView );
00066 
00067   // Copy the View to an vector for easier iterator manipulation convenience
00068   vector<const Candidate* > particles;
00069   for( CandidateView::const_iterator p = particlesViewH->begin();  p != particlesViewH->end(); ++p ) {
00070     particles.push_back(&*p);
00071   }
00072 
00073   // List of indices for the partons to add
00074   vector<int> partonIndices;
00075   // List of progenitors for those partons
00076   vector<int> progenitorIndices;
00077   // List of sisters for those partons
00078   vector<int> sisterIndices;
00079   // Flavor sources
00080   vector<FlavorHistory::FLAVOR_T> flavorSources;
00081   
00082   // Make a new flavor history vector
00083   auto_ptr<FlavorHistoryEvent > flavorHistoryEvent ( new FlavorHistoryEvent () ) ;
00084 
00085   // ------------------------------------------------------------
00086   // Loop over partons
00087   // ------------------------------------------------------------
00088   vector<const Candidate* >::size_type j;
00089   vector<const Candidate* >::size_type j_max = particles.size();
00090   for( j=0; j<j_max; ++j ) {
00091 
00092     // Get the candidate
00093     const Candidate *p = particles[j];
00094     // Set up indices that we'll need for the flavor history
00095     vector<Candidate const *>::size_type partonIndex=j;
00096     vector<Candidate const *>::size_type progenitorIndex=j_max;
00097     bool foundProgenitor = false; 
00098     FlavorHistory::FLAVOR_T flavorSource=FlavorHistory::FLAVOR_NULL;
00099 
00100 
00101     int idabs = std::abs( (p)->pdgId() );
00102     int nDa = (p)->numberOfDaughters();
00103 
00104     // Check if we have a status 2 or 3 particle, which is a parton before the string.
00105     // Only consider quarks. 
00106     if ( idabs <= FlavorHistory::tQuarkId && p->status() == 2 &&
00107          idabs == pdgIdToSelect_ ) {
00108       // Ensure the parton in question has daughters 
00109       if ( nDa > 0 ) {
00110         // Ensure the parton passes some minimum kinematic cuts
00111         if((p)->pt() > ptMinShower_ && fabs((p)->eta())<etaMaxShower_) {
00112 
00113           if(verbose_) cout << "--------------------------" << endl;
00114           if(verbose_) cout << "Processing particle " << j << ", status = " << p->status() << ", pdg id = " << p->pdgId() << endl;
00115 
00116 
00117           //Main (but simple) workhorse to get all ancestors
00118           vector<Candidate const *> allParents;
00119           getAncestors( *p, allParents );
00120             
00121           if(verbose_) {
00122             cout << "Parents : " << endl;
00123             vector<Candidate const *>::const_iterator iprint = allParents.begin(),
00124               iprintend = allParents.end();
00125             for ( ; iprint != iprintend; ++iprint ) 
00126               cout << " status = " << (*iprint)->status() << ", pdg id = " << (*iprint)->pdgId() << ", pt = " << (*iprint)->pt() << endl;
00127           }
00128           
00129           // boolean to check if there are status 3 partons of the same flavor.
00130           // If this is false, then the event is gluon splitting
00131           bool status3AncestorOfSameFlavor = false;
00132 
00133           // -------------------------------------------------------------
00134           // Identify the ancestry of the Quark
00135           // 
00136           // 
00137           // Matrix Element:
00138           //    Status 3 parent with precisely 2 "grandparents" that
00139           //    is outside of the "initial" section (0-5) that has the
00140           //    same ID as the status 2 parton in question. 
00141           //
00142           // Flavor excitation:
00143           //    If we find only one outgoing parton.
00144           //
00145           // Gluon splitting:
00146           //    Parent is a quark of a different flavor than the parton
00147           //    in question, or a gluon. 
00148           //    Can come from either ISR or FSR.
00149           //
00150           // True decay:
00151           //    Decays from a resonance like top, Higgs, etc.
00152           // -------------------------------------------------------------
00153           vector<Candidate const *>::size_type a_size = allParents.size();
00154 
00155           // 
00156           // Loop over all the ancestors of this parton and find the progenitor.
00157           // 
00158           for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
00159             const Candidate * aParent=allParents[i];
00160             vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
00161 
00162 
00163             // Get the index of the progenitor candidate
00164             progenitorIndex = found - particles.begin();
00165 
00166             int aParentId = std::abs(aParent->pdgId());
00167             
00168             // This will be used to check if there is gluon splitting present
00169             if ( aParent->status() == 3 && aParent->pdgId() == p->pdgId() ) status3AncestorOfSameFlavor = true;
00170 
00171             // -----------------------------------------------------------------------
00172             // Here we examine particles that were produced after the collision
00173             // -----------------------------------------------------------------------
00174             if( progenitorIndex > 5 ) {
00175               // Here is where we have a matrix element process.
00176               // The signature is to have a status 3 parent with precisely 2 parents
00177               // that is outside of the "initial" section (0-5) that has the
00178               // same ID as the status 2 parton in question.
00179               // NOTE: If this parton has no sister, then this will be classified as
00180               // a flavor excitation. Often the only difference is that the matrix element 
00181               // cases have a sister, while the flavor excitation cases do not. 
00182               // If we do not find a sister below, this will be classified as a flavor
00183               // excitation. 
00184               if(  aParent->numberOfMothers() == 2 &&  
00185                   aParent->pdgId() == p->pdgId() && aParent->status() == 3 ) {
00186                 if(verbose_) cout << "Matrix Element progenitor" << endl;
00187                 if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_ME;
00188                 foundProgenitor = true;
00189               } 
00190               // Here we have a true decay. Parent is not a quark or a gluon.
00191               else if( (aParentId>pdgIdToSelect_ && aParentId<FlavorHistory::gluonId) || 
00192                        aParentId > FlavorHistory::gluonId ) {
00193                 if(verbose_) cout << "Flavor decay progenitor" << endl;
00194                 if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_DECAY;
00195                 foundProgenitor = true;
00196               }
00197             }// end if progenitorIndex > 5
00198           }// end loop over parents
00199 
00200           // Otherwise, classify it as gluon splitting. Take the first status 3 parton with 1 parent
00201           // that you see as the progenitor
00202           if ( !foundProgenitor  ) {
00203             if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_GS;
00204             // Now get the parton with only one parent (the proton) and that is the progenitor
00205             for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
00206               const Candidate * aParent=allParents[i];
00207               vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
00208               // Get the index of the progenitor candidate
00209               progenitorIndex = found - particles.begin();
00210               
00211               if ( aParent->numberOfMothers() == 1 && aParent->status() == 3 ) {
00212                 foundProgenitor = true;
00213               }
00214             }// end loop over parents
00215           }// end if we haven't found a progenitor, and if we haven't found a status 3 parton of the same flavor
00216            // (i.e. end if examining gluon splitting)
00217         }// End if this parton passes some minimal kinematic cuts
00218       }// End if this particle is status 2 (has strings as daughters)
00219 
00220 
00221 
00222       // Make sure we've actually found a progenitor
00223       if ( !foundProgenitor ) progenitorIndex = j_max;
00224 
00225       // We've found the particle, add to the list
00226 
00227       partonIndices.push_back( partonIndex );
00228       progenitorIndices.push_back( progenitorIndex );
00229       flavorSources.push_back(flavorSource);
00230       sisterIndices.push_back( -1 ); // set below
00231         
00232     }// End if this particle was a status==2 parton
00233   }// end loop over particles
00234 
00235   // Now add sisters.
00236   // Also if the event is preliminarily classified as "matrix element", check to
00237   // make sure that they have a sister. If not, it is flavor excitation. 
00238   
00239   if ( verbose_ ) cout << "Making sisters" << endl;
00240   // First make sure nothing went terribly wrong:
00241   if ( partonIndices.size() == progenitorIndices.size() && partonIndices.size() > 0 ) {
00242     // Now loop over the candidates
00243     for ( unsigned int ii = 0; ii < partonIndices.size(); ++ii ) {
00244       // Get the iith particle
00245       const Candidate * candi = particles[partonIndices[ii]];
00246       if ( verbose_ ) cout << "   Sister candidate particle 1:  " << ii << ", pdgid = " << candi->pdgId() << ", status = " << candi->status() << endl;
00247       // Get the iith progenitor
00248       // Now loop over the other flavor history candidates and
00249       // attempt to find a sister to this one
00250       for ( unsigned int jj = 0; jj < partonIndices.size(); ++jj ) {
00251         if ( ii != jj ) {
00252           const Candidate * candj = particles[partonIndices[jj]];
00253       if ( verbose_ ) cout << "   Sister candidate particle 2:  " << jj << ", pdgid = " << candj->pdgId() << ", status = " << candj->status() << endl;
00254           // They should be opposite in pdgid and have the same status, and
00255           // the same progenitory.
00256           if ( candi->pdgId() == -1 * candj->pdgId() && candi->status() == candj->status()  ) {
00257             sisterIndices[ii] = partonIndices[jj];
00258             if ( verbose_ ) cout << "Parton " << partonIndices[ii] << " has sister " << sisterIndices[ii] << endl;
00259           }
00260         }
00261       }
00262 
00263       // Here, ensure that there is a sister. Otherwise this is "flavor excitation"
00264       if ( sisterIndices[ii] < 0 ) {
00265         if ( verbose_ ) cout << "No sister. Classified as flavor excitation" << endl;
00266         flavorSources[ii] = FlavorHistory::FLAVOR_EXC;
00267       } 
00268 
00269       if ( verbose_ ) cout << "Getting matched jet" << endl;
00270       // Get the closest match in the matched object collection
00271       CandidateView::const_iterator matched = getClosestJet( matchedView, *candi );
00272       CandidateView::const_iterator matchedBegin = matchedView->begin();
00273       CandidateView::const_iterator matchedEnd = matchedView->end();
00274 
00275       
00276       if ( verbose_ ) cout << "Getting sister jet" << endl;
00277       // Get the closest sister in the matched object collection, if sister is found
00278       CandidateView::const_iterator sister = 
00279         ( sisterIndices[ii] >= 0 && static_cast<unsigned int>(sisterIndices[ii]) < particles.size() ) ? 
00280         getClosestJet( matchedView, *particles[sisterIndices[ii]] ) :
00281         matchedEnd ;
00282 
00283       if ( verbose_ ) cout << "Making matched shallow clones : ";
00284       ShallowClonePtrCandidate matchedCand ;
00285       if ( matched != matchedEnd ) {
00286         if ( verbose_ ) cout << " found" << endl;
00287         matchedCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, matched - matchedBegin ) );
00288       } else {
00289         if ( verbose_ ) cout << " NOT found" << endl;
00290       }
00291       
00292       if ( verbose_ ) cout << "Making sister shallow clones : ";
00293       ShallowClonePtrCandidate sisterCand;
00294       if ( sister != matchedEnd ) {
00295         if ( verbose_ ) cout << " found" << endl;
00296         sisterCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, sister - matchedBegin ) );
00297       } else {
00298         if ( verbose_ ) cout << " NOT found" << endl;
00299       }
00300       
00301       if ( verbose_ ) cout << "Making history object" << endl;
00302       // Now create the flavor history object
00303       FlavorHistory history (flavorSources[ii], 
00304                              particlesViewH, 
00305                              partonIndices[ii], progenitorIndices[ii], sisterIndices[ii],
00306                              matchedCand,
00307                              sisterCand );
00308       if ( verbose_ ) cout << "Adding flavor history : " << history << endl;
00309       flavorHistoryEvent->push_back( history ); 
00310     }
00311   }
00312   
00313   // If we got any partons, cache them and then write them out
00314   if ( flavorHistoryEvent->size() > 0 ) {
00315 
00316     // Calculate some nice variables for FlavorHistoryEvent
00317     if ( verbose_ ) cout << "Caching flavor history event" << endl;
00318     flavorHistoryEvent->cache();
00319 
00320     if ( verbose_ ) {
00321       cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryEvent->size() << endl;
00322       vector<FlavorHistory>::const_iterator i = flavorHistoryEvent->begin(),
00323         iend = flavorHistoryEvent->end();
00324       for ( ; i !=iend; ++i ) {
00325         cout << *i << endl;
00326       }
00327     }
00328   }
00329 
00330   // Now add the flavor history to the event record
00331   evt.put( flavorHistoryEvent, flavorHistoryName_ );
00332 }
00333 
00334  
00335 // Helper function to get all ancestors of this candidate
00336 void FlavorHistoryProducer::getAncestors(const Candidate &c,
00337                                          vector<Candidate const *> & moms )
00338 {
00339 
00340   if( c.numberOfMothers() == 1 ) {
00341     const Candidate * dau = &c;
00342     const Candidate * mom = c.mother();
00343     while ( dau->numberOfMothers() != 0) {
00344       moms.push_back( dau );
00345       dau = mom ;
00346       mom = dau->mother();
00347     } 
00348   } 
00349 }
00350 
00351 
00352 CandidateView::const_iterator 
00353 FlavorHistoryProducer::getClosestJet( Handle<CandidateView> const & pJets,
00354                                       reco::Candidate const & parton ) const 
00355 {
00356   double dr = matchDR_;
00357   CandidateView::const_iterator j = pJets->begin(),
00358     jend = pJets->end();
00359   CandidateView::const_iterator bestJet = pJets->end();
00360   for ( ; j != jend; ++j ) {
00361     double dri = deltaR( parton.p4(), j->p4() );
00362     if ( dri < dr ) {
00363       dr = dri;
00364       bestJet = j;
00365     }
00366   }
00367   return bestJet;
00368 }
00369 
00370 #include "FWCore/Framework/interface/MakerMacros.h"
00371 
00372 DEFINE_FWK_MODULE( FlavorHistoryProducer );