JetMCUtils/src/JetMCTag.cc

0001 #include "PhysicsTools/JetMCUtils/interface/JetMCTag.h"
0002 #include "PhysicsTools/JetMCUtils/interface/CandMCTag.h"
0003 #include <iostream>
0004
0005 using namespace std;
0006 using namespace reco;
0007 using namespace JetMCTagUtils;
0008 using namespace CandMCTagUtils;
0009
0010 ///////////////////////////////////////////////////////////////////////
0011
0012 double JetMCTagUtils::EnergyRatioFromBHadrons(const Candidate &c) {
0013   double ratioForBjet = 0;
0014   double ratio = 0;
0015   for (Candidate::const_iterator itC = c.begin(); itC != c.end(); itC++) {
0016     bool isFromB = decayFromBHadron(*itC);
0017     ratio = itC->energy() / c.energy();
0018     if (isFromB)
0019       ratioForBjet += ratio;
0020   }
0021   return ratioForBjet;
0022 }
0023
0024 double JetMCTagUtils::EnergyRatioFromCHadrons(const Candidate &c) {
0025   double ratioForCjet = 0;
0026   double ratio = 0;
0027   for (Candidate::const_iterator itC = c.begin(); itC != c.end(); itC++) {
0028     bool isFromC = decayFromCHadron(*itC);
0029     ratio = itC->energy() / c.energy();
0030     if (isFromC)
0031       ratioForCjet += ratio;
0032   }
0033   return ratioForCjet;
0034 }
0035
0036 bool JetMCTagUtils::decayFromBHadron(const Candidate &c) {
0037   bool isFromB = false;
0038   vector<const Candidate *> allParents = getAncestors(c);
0039   for (vector<const Candidate *>::const_iterator aParent = allParents.begin(); aParent != allParents.end(); aParent++) {
0040     if (hasBottom(**aParent))
0041       isFromB = true;
0042     /*
0043          cout << "     particle Parent is " << (*aParent)->status()
0044               << " type " << (*aParent)->pdgId()
0045               << " pt=" << (*aParent)->pt()
0046               << " isB = " << isFromB
0047               << endl;
0048 */
0049   }
0050   return isFromB;
0051 }
0052
0053 bool JetMCTagUtils::decayFromCHadron(const Candidate &c) {
0054   bool isFromC = false;
0055   vector<const Candidate *> allParents = getAncestors(c);
0056   for (vector<const Candidate *>::const_iterator aParent = allParents.begin(); aParent != allParents.end(); aParent++) {
0057     if (hasCharm(**aParent))
0058       isFromC = true;
0059     /*
0060          cout << "     particle Parent is " << (*aParent)->status()
0061               << " type " << (*aParent)->pdgId()
0062               << " pt=" << (*aParent)->pt()
0063               << " isC = " << isFromC
0064               << endl;
0065 */
0066   }
0067   return isFromC;
0068 }
0069
0070 std::string JetMCTagUtils::genTauDecayMode(const CompositePtrCandidate &c) {
0071   int numElectrons = 0;
0072   int numMuons = 0;
0073   int numTaus = 0;
0074   int numChargedHadrons = 0;
0075   int numNeutralHadrons = 0;
0076   int numPhotons = 0;
0077
0078   const CompositePtrCandidate::daughters &daughters = c.daughterPtrVector();
0079   for (CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin(); daughter != daughters.end();
0080        ++daughter) {
0081     int pdg_id = abs((*daughter)->pdgId());
0082
0083     switch (pdg_id) {
0084       case 22:
0085         numPhotons++;
0086         break;
0087       case 11:
0088         numElectrons++;
0089         break;
0090       case 13:
0091         numMuons++;
0092         break;
0093       case 15:
0094         numTaus++;
0095         break;
0096       default: {
0097         if ((*daughter)->charge() != 0)
0098           numChargedHadrons++;
0099         else
0100           numNeutralHadrons++;
0101       }
0102     }
0103   }
0104
0105   if (numElectrons == 1)
0106     return std::string("electron");
0107   else if (numMuons == 1)
0108     return std::string("muon");
0109   else if (numTaus == 1)  //MB: a tau undecayed by generator or an intermediate state used to generate radiations
0110     return std::string("tau");
0111
0112   switch (numChargedHadrons) {
0113     case 1:
0114       if (numNeutralHadrons != 0)
0115         return std::string("oneProngOther");
0116       switch (numPhotons) {
0117         case 0:
0118           return std::string("oneProng0Pi0");
0119         case 2:
0120           return std::string("oneProng1Pi0");
0121         case 4:
0122           return std::string("oneProng2Pi0");
0123         default:
0124           return std::string("oneProngOther");
0125       }
0126     case 3:
0127       if (numNeutralHadrons != 0)
0128         return std::string("threeProngOther");
0129       switch (numPhotons) {
0130         case 0:
0131           return std::string("threeProng0Pi0");
0132         case 2:
0133           return std::string("threeProng1Pi0");
0134         default:
0135           return std::string("threeProngOther");
0136       }
0137     default:
0138       return std::string("rare");
0139   }
0140 }