ParticleFlow/src/FWPFTrackUtils.cc

0001 #include "Fireworks/ParticleFlow/interface/FWPFTrackUtils.h"
0002
0003 FWPFTrackSingleton *FWPFTrackSingleton::pInstance = nullptr;
0004 bool FWPFTrackSingleton::instanceFlag = false;
0005
0006 //______________________________________________________________________________
0007 FWPFTrackSingleton *FWPFTrackSingleton::Instance() {
0008   if (!instanceFlag)  // Instance doesn't exist yet
0009   {
0010     pInstance = new FWPFTrackSingleton();
0011     instanceFlag = true;
0012   }
0013
0014   return pInstance;  // Pointer to sole instance
0015 }
0016
0017 //______________________________________________________________________________
0018 void FWPFTrackSingleton::initPropagator() {
0019   m_magField = new FWMagField();
0020
0021   // Common propagator, helix stepper
0022   m_trackPropagator = new TEveTrackPropagator();
0023   m_trackPropagator->SetMagFieldObj(m_magField, false);
0024   m_trackPropagator->SetMaxR(FWPFGeom::caloR3());
0025   m_trackPropagator->SetMaxZ(FWPFGeom::caloZ2());
0026   m_trackPropagator->SetDelta(0.01);
0027   m_trackPropagator->SetProjTrackBreaking(TEveTrackPropagator::kPTB_UseLastPointPos);
0028   m_trackPropagator->SetRnrPTBMarkers(kTRUE);
0029   m_trackPropagator->IncDenyDestroy();
0030
0031   // Tracker propagator
0032   m_trackerTrackPropagator = new TEveTrackPropagator();
0033   m_trackerTrackPropagator->SetStepper(TEveTrackPropagator::kRungeKutta);
0034   m_trackerTrackPropagator->SetMagFieldObj(m_magField, false);
0035   m_trackerTrackPropagator->SetDelta(0.01);
0036   m_trackerTrackPropagator->SetMaxR(FWPFGeom::caloR3());
0037   m_trackerTrackPropagator->SetMaxZ(FWPFGeom::caloZ2());
0038   m_trackerTrackPropagator->SetProjTrackBreaking(TEveTrackPropagator::kPTB_UseLastPointPos);
0039   m_trackerTrackPropagator->SetRnrPTBMarkers(kTRUE);
0040   m_trackerTrackPropagator->IncDenyDestroy();
0041 }
0042
0043 //______________________________________________________________________________
0044 FWPFTrackUtils::FWPFTrackUtils() { m_singleton = FWPFTrackSingleton::Instance(); }
0045
0046 //______________________________________________________________________________
0047 TEveTrack *FWPFTrackUtils::getTrack(const reco::Track &iData) {
0048   TEveTrackPropagator *propagator =
0049       (!iData.extra().isAvailable()) ? m_singleton->getTrackerTrackPropagator() : m_singleton->getTrackPropagator();
0050
0051   TEveRecTrack t;
0052   t.fBeta = 1;
0053   t.fP = TEveVector(iData.px(), iData.py(), iData.pz());
0054   t.fV = TEveVector(iData.vertex().x(), iData.vertex().y(), iData.vertex().z());
0055   t.fSign = iData.charge();
0056   TEveTrack *trk = new TEveTrack(&t, propagator);
0057   trk->MakeTrack();
0058
0059   return trk;
0060 }
0061
0062 //______________________________________________________________________________
0063 TEveStraightLineSet *FWPFTrackUtils::setupLegoTrack(const reco::Track &iData) {
0064   using namespace FWPFGeom;
0065
0066   // Declarations
0067   int wraps[3] = {-1, -1, -1};
0068   bool ECAL = false;
0069   TEveTrack *trk = getTrack(iData);
0070   std::vector<TEveVector> trackPoints(trk->GetN() - 1);
0071   const Float_t *points = trk->GetP();
0072   TEveStraightLineSet *legoTrack = new TEveStraightLineSet();
0073
0074   if (m_singleton->getField()->getSource() == FWMagField::kNone) {
0075     if (fabs(iData.eta()) < 2.0 && iData.pt() > 0.5 && iData.pt() < 30) {
0076       double estimate = fw::estimate_field(iData, true);
0077       if (estimate >= 0)
0078         m_singleton->getField()->guessField(estimate);
0079     }
0080   }
0081
0082   // Convert to Eta/Phi and store in vector
0083   for (Int_t i = 1; i < trk->GetN(); ++i) {
0084     int j = i * 3;
0085     TEveVector temp = TEveVector(points[j], points[j + 1], points[j + 2]);
0086     TEveVector vec = TEveVector(temp.Eta(), temp.Phi(), 0.001);
0087
0088     trackPoints[i - 1] = vec;
0089   }
0090
0091   // Add first point to ps if necessary
0092   for (Int_t i = 1; i < trk->GetN(); ++i) {
0093     int j = i * 3;
0094     TEveVector v1 = TEveVector(points[j], points[j + 1], points[j + 2]);
0095
0096     if (!ECAL) {
0097       if (FWPFMaths::checkIntersect(v1, caloR1())) {
0098         TEveVector v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0099         TEveVector xyPoint = FWPFMaths::lineCircleIntersect(v1, v2, caloR1());
0100         TEveVector zPoint;
0101         if (v1.fZ < 0)
0102           zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ - 50.f);
0103         else
0104           zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ + 50.f);
0105
0106         TEveVector vec = FWPFMaths::lineLineIntersect(v1, v2, xyPoint, zPoint);
0107         legoTrack->AddMarker(vec.Eta(), vec.Phi(), 0.001, 0);
0108
0109         wraps[0] = i;  // There is now a chance that the track will also reach the HCAL radius
0110         ECAL = true;
0111       } else if (fabs(v1.fZ) >= caloZ1()) {
0112         TEveVector p1, p2;
0113         TEveVector vec, v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0114         float z, y = FWPFMaths::linearInterpolation(v2, v1, caloZ1());
0115
0116         if (v2.fZ > 0)
0117           z = caloZ1();
0118         else
0119           z = caloZ1() * -1;
0120
0121         p1 = TEveVector(v2.fX + 50, y, z);
0122         p2 = TEveVector(v2.fX - 50, y, z);
0123         vec = FWPFMaths::lineLineIntersect(v1, v2, p1, p2);
0124
0125         legoTrack->AddMarker(vec.Eta(), vec.Phi(), 0.001, 0);
0126         wraps[0] = i;
0127         ECAL = true;
0128       }
0129     } else if (FWPFMaths::checkIntersect(v1, caloR2())) {
0130       TEveVector v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0131       TEveVector xyPoint = FWPFMaths::lineCircleIntersect(v1, v2, caloR2());
0132       TEveVector zPoint;
0133       if (v1.fZ < 0)
0134         zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ - 50.f);
0135       else
0136         zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ + 50.f);
0137
0138       TEveVector vec = FWPFMaths::lineLineIntersect(v1, v2, xyPoint, zPoint);
0139       legoTrack->AddMarker(vec.Eta(), vec.Phi(), 0.001, 1);
0140
0141       wraps[1] = i;  // There is now a chance that the track will also reach the HCAL radius
0142       break;
0143     }
0144   }
0145
0146   if (wraps[0] != -1)  //if( ECAL )
0147   {
0148     int i = (trk->GetN() - 1) * 3;
0149     int j = trk->GetN() - 2;  // This is equal to the last index in trackPoints vector
0150     TEveVector vec = TEveVector(points[i], points[i + 1], points[i + 2]);
0151
0152     if (FWPFMaths::checkIntersect(vec, caloR3() - 1)) {
0153       legoTrack->AddMarker(vec.Eta(), vec.Phi(), 0.001, 2);
0154       wraps[2] = j;
0155     } else if (fabs(vec.fZ) >= caloZ2()) {
0156       legoTrack->AddMarker(vec.Eta(), vec.Phi(), 0.001, 2);
0157       wraps[2] = j;
0158     }
0159   }
0160
0161   /* Handle phi wrapping */
0162   for (int i = 0; i < static_cast<int>(trackPoints.size() - 1); ++i) {
0163     if ((trackPoints[i + 1].fY - trackPoints[i].fY) > 1) {
0164       trackPoints[i + 1].fY -= TMath::TwoPi();
0165       if (i == wraps[0]) {
0166         TEveChunkManager::iterator mi(legoTrack->GetMarkerPlex());
0167         mi.next();  // First point
0168         TEveStraightLineSet::Marker_t &m = *(TEveStraightLineSet::Marker_t *)mi();
0169         m.fV[0] = trackPoints[i + 1].fX;
0170         m.fV[1] = trackPoints[i + 1].fY;
0171         m.fV[2] = 0.001;
0172       } else if (i == wraps[1]) {
0173         TEveChunkManager::iterator mi(legoTrack->GetMarkerPlex());
0174         mi.next();
0175         mi.next();  // Second point
0176         TEveStraightLineSet::Marker_t &m = *(TEveStraightLineSet::Marker_t *)mi();
0177         m.fV[0] = trackPoints[i + 1].fX;
0178         m.fV[1] = trackPoints[i + 1].fY;
0179         m.fV[2] = 0.001;
0180       }
0181     }
0182
0183     if ((trackPoints[i].fY - trackPoints[i + 1].fY) > 1) {
0184       trackPoints[i + 1].fY += TMath::TwoPi();
0185       if (i == wraps[0]) {
0186         TEveChunkManager::iterator mi(legoTrack->GetMarkerPlex());
0187         mi.next();  // First point
0188         TEveStraightLineSet::Marker_t &m = *(TEveStraightLineSet::Marker_t *)mi();
0189         m.fV[0] = trackPoints[i + 1].fX;
0190         m.fV[1] = trackPoints[i + 1].fY;
0191         m.fV[2] = 0.001;
0192       } else if (i == wraps[1]) {
0193         TEveChunkManager::iterator mi(legoTrack->GetMarkerPlex());
0194         mi.next();
0195         mi.next();  // Second point
0196         TEveStraightLineSet::Marker_t &m = *(TEveStraightLineSet::Marker_t *)mi();
0197         m.fV[0] = trackPoints[i + 1].fX;
0198         m.fV[1] = trackPoints[i + 1].fY;
0199         m.fV[2] = 0.001;
0200       }
0201     }
0202   }
0203
0204   int end = static_cast<int>(trackPoints.size() - 1);
0205   if (wraps[2] == end) {
0206     TEveChunkManager::iterator mi(legoTrack->GetMarkerPlex());
0207     mi.next();
0208     mi.next();
0209     mi.next();  // Third point
0210     TEveStraightLineSet::Marker_t &m = *(TEveStraightLineSet::Marker_t *)mi();
0211     m.fV[0] = trackPoints[end].fX;
0212     m.fV[1] = trackPoints[end].fY;
0213     m.fV[2] = 0.001;
0214   }
0215
0216   // Set points on TEveLineSet object ready for displaying
0217   for (unsigned int i = 0; i < trackPoints.size() - 1; ++i)
0218     legoTrack->AddLine(trackPoints[i], trackPoints[i + 1]);
0219
0220   legoTrack->SetDepthTest(false);
0221   legoTrack->SetMarkerStyle(4);
0222   legoTrack->SetMarkerSize(1);
0223   legoTrack->SetRnrMarkers(true);
0224
0225   delete trk;  // Release memory that is no longer required
0226
0227   return legoTrack;
0228 }
0229
0230 //______________________________________________________________________________
0231 TEveTrack *FWPFTrackUtils::setupTrack(const reco::Track &iData) {
0232   if (m_singleton->getField()->getSource() == FWMagField::kNone) {
0233     if (fabs(iData.eta()) < 2.0 && iData.pt() > 0.5 && iData.pt() < 30) {
0234       double estimate = fw::estimate_field(iData, true);
0235       if (estimate >= 0)
0236         m_singleton->getField()->guessField(estimate);
0237     }
0238   }
0239
0240   TEveTrack *trk = getTrack(iData);
0241
0242   return trk;
0243 }
0244
0245 //______________________________________________________________________________
0246 TEvePointSet *FWPFTrackUtils::getCollisionMarkers(const TEveTrack *trk) {
0247   using namespace FWPFGeom;
0248
0249   bool ECAL = false;
0250   const Float_t *points = trk->GetP();
0251   TEvePointSet *ps = new TEvePointSet();
0252
0253   for (Int_t i = 1; i < trk->GetN(); ++i) {
0254     int j = i * 3;
0255     TEveVector v1 = TEveVector(points[j], points[j + 1], points[j + 2]);
0256
0257     if (!ECAL) {
0258       if (FWPFMaths::checkIntersect(v1, caloR1())) {
0259         TEveVector v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0260         TEveVector xyPoint = FWPFMaths::lineCircleIntersect(v1, v2, caloR1());
0261         TEveVector zPoint;
0262         if (v1.fZ < 0)
0263           zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ - 50.f);
0264         else
0265           zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ + 50.f);
0266
0267         TEveVector vec = FWPFMaths::lineLineIntersect(v1, v2, xyPoint, zPoint);
0268         ps->SetNextPoint(vec.fX, vec.fY, vec.fZ);
0269
0270         ECAL = true;
0271       } else if (fabs(v1.fZ) >= caloZ1()) {
0272         TEveVector p1, p2;
0273         TEveVector vec, v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0274         float z, y = FWPFMaths::linearInterpolation(v2, v1, caloZ1());
0275
0276         if (v2.fZ > 0)
0277           z = caloZ1();
0278         else
0279           z = caloZ1() * -1;
0280
0281         p1 = TEveVector(v2.fX + 50, y, z);
0282         p2 = TEveVector(v2.fX - 50, y, z);
0283         vec = FWPFMaths::lineLineIntersect(v1, v2, p1, p2);
0284
0285         ps->SetNextPoint(vec.fX, vec.fY, vec.fZ);
0286         ECAL = true;
0287       }
0288     } else if (FWPFMaths::checkIntersect(v1, caloR2())) {
0289       TEveVector v2 = TEveVector(points[j - 3], points[j - 2], points[j - 1]);
0290       TEveVector xyPoint = FWPFMaths::lineCircleIntersect(v1, v2, caloR2());
0291       TEveVector zPoint;
0292       if (v1.fZ < 0)
0293         zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ - 50.f);
0294       else
0295         zPoint = TEveVector(xyPoint.fX, xyPoint.fY, v1.fZ + 50.f);
0296
0297       TEveVector vec = FWPFMaths::lineLineIntersect(v1, v2, xyPoint, zPoint);
0298       ps->SetNextPoint(vec.fX, vec.fY, vec.fZ);
0299       break;  // ECAL and HCAL collisions found so stop looping
0300     }
0301   }
0302
0303   // HCAL collisions (outer radius and endcap)
0304   int i = (trk->GetN() - 1) * 3;
0305   TEveVector vec = TEveVector(points[i], points[i + 1], points[i + 2]);
0306
0307   if (FWPFMaths::checkIntersect(vec, caloR3() - 1))
0308     ps->SetNextPoint(vec.fX, vec.fY, vec.fZ);
0309   else if (fabs(vec.fZ) >= caloZ2())
0310     ps->SetNextPoint(vec.fX, vec.fY, vec.fZ);
0311
0312   return ps;
0313 }