ParticleFlowReco/src/PFDisplacedVertex.cc

0001 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
0002
0003 #include "TMath.h"
0004
0005 using namespace std;
0006 using namespace reco;
0007
0008 PFDisplacedVertex::PFDisplacedVertex() : Vertex(), vertexType_(ANY), primaryDirection_(0, 0, 0) {}
0009
0010 PFDisplacedVertex::PFDisplacedVertex(Vertex& v) : Vertex(v), vertexType_(ANY), primaryDirection_(0, 0, 0) {}
0011
0012 void PFDisplacedVertex::addElement(const TrackBaseRef& r,
0013                                    const Track& refTrack,
0014                                    const PFTrackHitFullInfo& hitInfo,
0015                                    VertexTrackType trackType,
0016                                    float w) {
0017   add(r, refTrack, w);
0018   trackTypes_.push_back(trackType);
0019   trackHitFullInfos_.push_back(hitInfo);
0020 }
0021
0022 void PFDisplacedVertex::cleanTracks() {
0023   removeTracks();
0024   trackTypes_.clear();
0025   trackHitFullInfos_.clear();
0026 }
0027
0028 const bool PFDisplacedVertex::isThereKindTracks(VertexTrackType T) const {
0029   vector<VertexTrackType>::const_iterator iter = find(trackTypes_.begin(), trackTypes_.end(), T);
0030   return (iter != trackTypes_.end());
0031 }
0032
0033 const int PFDisplacedVertex::nKindTracks(VertexTrackType T) const {
0034   return count(trackTypes_.begin(), trackTypes_.end(), T);
0035 }
0036
0037 const size_t PFDisplacedVertex::trackPosition(const reco::TrackBaseRef& originalTrack) const {
0038   size_t pos = -1;
0039
0040   const Track refittedTrack = PFDisplacedVertex::refittedTrack(originalTrack);
0041
0042   std::vector<Track> refitTrks = refittedTracks();
0043   for (size_t i = 0; i < refitTrks.size(); i++) {
0044     if (fabs(refitTrks[i].pt() - refittedTrack.pt()) < 1.e-5) {
0045       pos = i;
0046       continue;
0047     }
0048   }
0049   //  cout << "pos = " << pos << endl;
0050
0051   return pos;
0052 }
0053
0054 void PFDisplacedVertex::setPrimaryDirection(const math::XYZPoint& pvtx) {
0055   primaryDirection_ = math::XYZVector(position().x(), position().y(), position().z());
0056   math::XYZVector vtx(pvtx.x(), pvtx.y(), pvtx.z());
0057
0058   primaryDirection_ = primaryDirection_ - vtx;
0059   primaryDirection_ /= (sqrt(primaryDirection_.Mag2()) + 1e-10);
0060 }
0061
0062 std::string PFDisplacedVertex::nameVertexType() const {
0063   switch (vertexType_) {
0064     case ANY:
0065       return "ANY";
0066     case FAKE:
0067       return "FAKE";
0068     case LOOPER:
0069       return "LOOPER";
0070     case NUCL:
0071       return "NUCL";
0072     case NUCL_LOOSE:
0073       return "NUCL_LOOSE";
0074     case NUCL_KINK:
0075       return "NUCL_KINK";
0076     case CONVERSION:
0077       return "CONVERSION";
0078     case CONVERSION_LOOSE:
0079       return "CONVERSION_LOOSE";
0080     case CONVERTED_BREMM:
0081       return "CONVERTED_BREMM";
0082     case K0_DECAY:
0083       return "K0_DECAY";
0084     case LAMBDA_DECAY:
0085       return "LAMBDA_DECAY";
0086     case LAMBDABAR_DECAY:
0087       return "LAMBDABAR_DECAY";
0088     case KPLUS_DECAY:
0089       return "KPLUS_DECAY";
0090     case KMINUS_DECAY:
0091       return "KMINUS_DECAY";
0092     case KPLUS_DECAY_LOOSE:
0093       return "KPLUS_DECAY_LOOSE";
0094     case KMINUS_DECAY_LOOSE:
0095       return "KMINUS_DECAY_LOOSE";
0096     case BSM_VERTEX:
0097       return "BSM_VERTEX";
0098     default:
0099       return "?";
0100   }
0101   return "?";
0102 }
0103
0104 const math::XYZTLorentzVector PFDisplacedVertex::momentum(string massHypo,
0105                                                           VertexTrackType T,
0106                                                           bool useRefitted,
0107                                                           double mass) const {
0108   M_Hypo mHypo = M_CUSTOM;
0109
0110   if (massHypo.find("PI") != string::npos)
0111     mHypo = M_PION;
0112   else if (massHypo.find("KAON") != string::npos)
0113     mHypo = M_KAON;
0114   else if (massHypo.find("LAMBDA") != string::npos)
0115     mHypo = M_LAMBDA;
0116   else if (massHypo.find("MASSLESS") != string::npos)
0117     mHypo = M_MASSLESS;
0118   else if (massHypo.find("CUSTOM") != string::npos)
0119     mHypo = M_CUSTOM;
0120
0121   return momentum(mHypo, T, useRefitted, mass);
0122 }
0123
0124 const math::XYZTLorentzVector PFDisplacedVertex::momentum(M_Hypo massHypo,
0125                                                           VertexTrackType T,
0126                                                           bool useRefitted,
0127                                                           double mass) const {
0128   const double m2 = getMass2(massHypo, mass);
0129
0130   math::XYZTLorentzVector P;
0131
0132   for (size_t i = 0; i < tracksSize(); i++) {
0133     bool bType = (trackTypes_[i] == T);
0134     if (T == T_TO_VERTEX || T == T_MERGED)
0135       bType = (trackTypes_[i] == T_TO_VERTEX || trackTypes_[i] == T_MERGED);
0136
0137     if (bType) {
0138       if (!useRefitted) {
0139         TrackBaseRef trackRef = originalTrack(refittedTracks()[i]);
0140
0141         double p2 = trackRef->momentum().Mag2();
0142         P += math::XYZTLorentzVector(
0143             trackRef->momentum().x(), trackRef->momentum().y(), trackRef->momentum().z(), sqrt(m2 + p2));
0144       } else {
0145         //  cout << "m2 " << m2 << endl;
0146
0147         double p2 = refittedTracks()[i].momentum().Mag2();
0148         P += math::XYZTLorentzVector(refittedTracks()[i].momentum().x(),
0149                                      refittedTracks()[i].momentum().y(),
0150                                      refittedTracks()[i].momentum().z(),
0151                                      sqrt(m2 + p2));
0152       }
0153     }
0154   }
0155
0156   return P;
0157 }
0158
0159 const int PFDisplacedVertex::totalCharge() const {
0160   int charge = 0;
0161
0162   for (size_t i = 0; i < tracksSize(); i++) {
0163     if (trackTypes_[i] == T_TO_VERTEX)
0164       charge += refittedTracks()[i].charge();
0165     else if (trackTypes_[i] == T_FROM_VERTEX)
0166       charge -= refittedTracks()[i].charge();
0167   }
0168
0169   return charge;
0170 }
0171
0172 const double PFDisplacedVertex::angle_io() const {
0173   math::XYZTLorentzVector momentumSec = secondaryMomentum((string) "PI", true);
0174
0175   math::XYZVector p_out = momentumSec.Vect();
0176
0177   math::XYZVector p_in = primaryDirection();
0178
0179   if (p_in.Mag2() < 1e-10)
0180     return -1;
0181   return acos(p_in.Dot(p_out) / sqrt(p_in.Mag2() * p_out.Mag2())) / TMath::Pi() * 180.0;
0182 }
0183
0184 const math::XYZVector PFDisplacedVertex::primaryDirection() const {
0185   math::XYZTLorentzVector momentumPrim = primaryMomentum((string) "PI", true);
0186   math::XYZTLorentzVector momentumSec = secondaryMomentum((string) "PI", true);
0187
0188   math::XYZVector p_in;
0189
0190   if ((isThereKindTracks(T_TO_VERTEX) || isThereKindTracks(T_MERGED)) && momentumPrim.E() > momentumSec.E()) {
0191     p_in = momentumPrim.Vect() / sqrt(momentumPrim.Vect().Mag2() + 1e-10);
0192   } else {
0193     p_in = primaryDirection_;
0194   }
0195
0196   return p_in;
0197 }
0198
0199 const double PFDisplacedVertex::getMass2(M_Hypo massHypo, double mass) const {
0200   // pion_mass = 0.1396 GeV
0201   double pion_mass2 = 0.0194;
0202   // k0_mass = 0.4976 GeV
0203   double kaon_mass2 = 0.2476;
0204   // lambda0_mass = 1.116 GeV
0205   double lambda_mass2 = 1.267;
0206
0207   if (massHypo == M_PION)
0208     return pion_mass2;
0209   else if (massHypo == M_KAON)
0210     return kaon_mass2;
0211   else if (massHypo == M_LAMBDA)
0212     return lambda_mass2;
0213   else if (massHypo == M_MASSLESS)
0214     return 0;
0215   else if (massHypo == M_CUSTOM)
0216     return mass * mass;
0217
0218   cout << "Warning: undefined mass hypothesis" << endl;
0219   return 0;
0220 }
0221
0222 void PFDisplacedVertex::Dump(ostream& out) const {
0223   if (!out)
0224     return;
0225
0226   out << "" << endl;
0227   out << "==================== This is a Displaced Vertex type " << nameVertexType() << " ===============" << endl;
0228
0229   out << " Vertex chi2 = " << chi2() << " ndf = " << ndof() << " normalised chi2 = " << normalizedChi2() << endl;
0230
0231   out << " The vertex Fitted Position is: x = " << position().x() << " y = " << position().y()
0232       << " rho = " << position().rho() << " z = " << position().z() << endl;
0233
0234   out << "\t--- Structure ---  " << endl;
0235   out << "Number of tracks: " << nTracks() << " nPrimary " << nPrimaryTracks() << " nMerged " << nMergedTracks()
0236       << " nSecondary " << nSecondaryTracks() << endl;
0237
0238   vector<PFDisplacedVertex::PFTrackHitFullInfo> pattern = trackHitFullInfos();
0239   vector<PFDisplacedVertex::VertexTrackType> trackType = trackTypes();
0240   for (unsigned i = 0; i < pattern.size(); i++) {
0241     out << "track " << i << " type = " << trackType[i] << " nHit BeforeVtx = " << pattern[i].first.first
0242         << " AfterVtx = " << pattern[i].second.first << " MissHit BeforeVtx = " << pattern[i].first.second
0243         << " AfterVtx = " << pattern[i].second.second << endl;
0244   }
0245
0246   math::XYZTLorentzVector mom_prim = primaryMomentum((string) "PI", true);
0247   math::XYZTLorentzVector mom_sec = secondaryMomentum((string) "PI", true);
0248
0249   // out << "Primary P:\t E " << setprecision(3) << setw(5) << mom_prim.E()
0250   out << "Primary P:\t E " << mom_prim.E() << "\tPt = " << mom_prim.Pt() << "\tPz = " << mom_prim.Pz()
0251       << "\tM = " << mom_prim.M() << "\tEta = " << mom_prim.Eta() << "\tPhi = " << mom_prim.Phi() << endl;
0252
0253   out << "Secondary P:\t E " << mom_sec.E() << "\tPt = " << mom_sec.Pt() << "\tPz = " << mom_sec.Pz()
0254       << "\tM = " << mom_sec.M() << "\tEta = " << mom_sec.Eta() << "\tPhi = " << mom_sec.Phi() << endl;
0255
0256   out << " The vertex Direction is x = " << primaryDirection().x() << " y = " << primaryDirection().y()
0257       << " z = " << primaryDirection().z() << " eta = " << primaryDirection().eta()
0258       << " phi = " << primaryDirection().phi() << endl;
0259
0260   out << " Angle_io = " << angle_io() << " deg" << endl << endl;
0261 }