Muons/src/FWMuonBuilder.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     Muons
0004 // Class  :     FWMuonBuilder
0005 //
0006
0007 #include "TEveVSDStructs.h"
0008 #include "TEveTrack.h"
0009 #include "TEveStraightLineSet.h"
0010 #include "TEveGeoNode.h"
0011 #include "TEveGeoShape.h"
0012 #include "TGeoArb8.h"
0013
0014 #include "Fireworks/Core/interface/FWEventItem.h"
0015 #include "Fireworks/Core/interface/FWMagField.h"
0016 #include "Fireworks/Core/interface/FWProxyBuilderBase.h"
0017 #include "Fireworks/Core/interface/FWGeometry.h"
0018 #include "Fireworks/Core/interface/fwLog.h"
0019
0020 #include "Fireworks/Candidates/interface/CandidateUtils.h"
0021
0022 #include "Fireworks/Tracks/interface/TrackUtils.h"
0023 #include "Fireworks/Tracks/interface/estimate_field.h"
0024
0025 #include "Fireworks/Muons/interface/FWMuonBuilder.h"
0026 #include "Fireworks/Muons/interface/SegmentUtils.h"
0027
0028 #include "DataFormats/MuonReco/interface/Muon.h"
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0031
0032 namespace {
0033   std::vector<TEveVector> getRecoTrajectoryPoints(const reco::Muon* muon, const FWEventItem* iItem) {
0034     std::vector<TEveVector> points;
0035     const FWGeometry* geom = iItem->getGeom();
0036
0037     float localTrajectoryPoint[3];
0038     float globalTrajectoryPoint[3];
0039
0040     const std::vector<reco::MuonChamberMatch>& matches = muon->matches();
0041     for (std::vector<reco::MuonChamberMatch>::const_iterator chamber = matches.begin(), chamberEnd = matches.end();
0042          chamber != chamberEnd;
0043          ++chamber) {
0044       // expected track position
0045       localTrajectoryPoint[0] = chamber->x;
0046       localTrajectoryPoint[1] = chamber->y;
0047       localTrajectoryPoint[2] = 0;
0048
0049       unsigned int rawid = chamber->id.rawId();
0050       if (geom->contains(rawid)) {
0051         geom->localToGlobal(rawid, localTrajectoryPoint, globalTrajectoryPoint);
0052         points.push_back(TEveVector(globalTrajectoryPoint[0], globalTrajectoryPoint[1], globalTrajectoryPoint[2]));
0053       }
0054     }
0055     return points;
0056   }
0057
0058   //______________________________________________________________________________
0059
0060   void addMatchInformation(const reco::Muon* muon, FWProxyBuilderBase* pb, TEveElement* parentList, bool showEndcap) {
0061     std::set<unsigned int> ids;
0062     const FWGeometry* geom = pb->context().getGeom();
0063
0064     const std::vector<reco::MuonChamberMatch>& matches = muon->matches();
0065
0066     //need to use auto_ptr since the segmentSet may not be passed to muonList
0067     std::unique_ptr<TEveStraightLineSet> segmentSet(new TEveStraightLineSet);
0068     // FIXME: This should be set elsewhere.
0069     segmentSet->SetLineWidth(4);
0070
0071     for (std::vector<reco::MuonChamberMatch>::const_iterator chamber = matches.begin(), chambersEnd = matches.end();
0072          chamber != chambersEnd;
0073          ++chamber) {
0074       unsigned int rawid = chamber->id.rawId();
0075       float segmentLength = 0.0;
0076       float segmentLimit = 0.0;
0077
0078       if (geom->contains(rawid)) {
0079         TEveGeoShape* shape = geom->getEveShape(rawid);
0080         shape->SetElementName("Chamber");
0081         shape->RefMainTrans().Scale(0.999, 0.999, 0.999);
0082
0083         FWGeometry::IdToInfoItr det = geom->find(rawid);
0084         if (det->shape[0] == 1)  // TGeoTrap
0085         {
0086           segmentLength = det->shape[3];
0087           segmentLimit = det->shape[4];
0088         } else if (det->shape[0] == 2)  // TGeoBBox
0089         {
0090           segmentLength = det->shape[3];
0091         } else {
0092           const double segmentLength = 15;
0093           fwLog(fwlog::kWarning) << Form(
0094               "FWMuonBuilder: unknown shape type in muon chamber with detId=%d. Setting segment length to %.0f cm.\n",
0095               rawid,
0096               segmentLength);
0097         }
0098
0099         if (ids.insert(rawid).second &&  // ensure that we add same chamber only once
0100             (chamber->detector() != MuonSubdetId::CSC || showEndcap)) {
0101           pb->setupAddElement(shape, parentList);
0102         }
0103
0104         for (std::vector<reco::MuonSegmentMatch>::const_iterator segment = chamber->segmentMatches.begin(),
0105                                                                  segmentEnd = chamber->segmentMatches.end();
0106              segment != segmentEnd;
0107              ++segment) {
0108           float segmentPosition[3] = {segment->x, segment->y, 0.0};
0109           float segmentDirection[3] = {segment->dXdZ, segment->dYdZ, 0.0};
0110
0111           float localSegmentInnerPoint[3];
0112           float localSegmentOuterPoint[3];
0113
0114           fireworks::createSegment(chamber->detector(),
0115                                    true,
0116                                    segmentLength,
0117                                    segmentLimit,
0118                                    segmentPosition,
0119                                    segmentDirection,
0120                                    localSegmentInnerPoint,
0121                                    localSegmentOuterPoint);
0122
0123           float globalSegmentInnerPoint[3];
0124           float globalSegmentOuterPoint[3];
0125
0126           geom->localToGlobal(*det, localSegmentInnerPoint, globalSegmentInnerPoint);
0127           geom->localToGlobal(*det, localSegmentOuterPoint, globalSegmentOuterPoint);
0128
0129           segmentSet->AddLine(globalSegmentInnerPoint[0],
0130                               globalSegmentInnerPoint[1],
0131                               globalSegmentInnerPoint[2],
0132                               globalSegmentOuterPoint[0],
0133                               globalSegmentOuterPoint[1],
0134                               globalSegmentOuterPoint[2]);
0135         }
0136       }
0137     }
0138
0139     if (!matches.empty())
0140       pb->setupAddElement(segmentSet.release(), parentList);
0141   }
0142
0143   //______________________________________________________________________________
0144
0145   bool buggyMuon(const reco::Muon* muon, const FWGeometry* geom) {
0146     if (!muon->standAloneMuon().isAvailable() || !muon->standAloneMuon()->extra().isAvailable())
0147       return false;
0148
0149     float localTrajectoryPoint[3];
0150     float globalTrajectoryPoint[3];
0151
0152     const std::vector<reco::MuonChamberMatch>& matches = muon->matches();
0153     for (std::vector<reco::MuonChamberMatch>::const_iterator chamber = matches.begin(), chamberEnd = matches.end();
0154          chamber != chamberEnd;
0155          ++chamber) {
0156       // expected track position
0157       localTrajectoryPoint[0] = chamber->x;
0158       localTrajectoryPoint[1] = chamber->y;
0159       localTrajectoryPoint[2] = 0;
0160
0161       unsigned int rawid = chamber->id.rawId();
0162       if (geom->contains(rawid)) {
0163         geom->localToGlobal(rawid, localTrajectoryPoint, globalTrajectoryPoint);
0164         double phi = atan2(globalTrajectoryPoint[1], globalTrajectoryPoint[0]);
0165         if (cos(phi - muon->standAloneMuon()->innerPosition().phi()) < 0)
0166           return true;
0167       }
0168     }
0169     return false;
0170   }
0171
0172   TEveTrack* prepareMuonTrackWithExtraPoints(const reco::Track& track,
0173                                              TEveTrackPropagator* propagator,
0174                                              const std::vector<TEveVector>& extraPoints) {
0175     TEveRecTrack t;
0176     t.fBeta = 1.;
0177     t.fSign = track.charge();
0178     t.fV.Set(track.vx(), track.vy(), track.vz());
0179     t.fP.Set(track.px(), track.py(), track.pz());
0180     //  t.fSign = muon->charge();
0181     //  t.fV.Set(muon->vx(), muon->vy(), muon->vz());
0182     //  t.fP.Set(muon->px(), muon->py(), muon->pz());
0183     TEveTrack* trk = new TEveTrack(&t, propagator);
0184     size_t n = extraPoints.size();
0185
0186     if (n > 1) {
0187       int lastDaughter = n - 2;
0188       for (int i = 0; i <= lastDaughter; ++i)
0189         trk->AddPathMark(TEvePathMark(TEvePathMark::kDaughter, extraPoints[i]));
0190     }
0191     trk->AddPathMark(TEvePathMark(TEvePathMark::kDecay, extraPoints.back()));
0192     return trk;
0193   }
0194
0195 }  // namespace
0196
0197 //
0198 // constructors and destructor
0199 //
0200 FWMuonBuilder::FWMuonBuilder() : m_lineWidth(1) {}
0201
0202 FWMuonBuilder::~FWMuonBuilder() {}
0203
0204 //
0205 // member functions
0206 //
0207 //______________________________________________________________________________
0208
0209 void FWMuonBuilder::calculateField(const reco::Muon& iData, FWMagField* field) {
0210   // if auto field estimation mode, do extra loop over muons.
0211   // use both inner and outer track if available
0212   if (field->getSource() == FWMagField::kNone) {
0213     if (fabs(iData.eta()) > 2.0 || iData.pt() < 3)
0214       return;
0215     if (iData.innerTrack().isAvailable()) {
0216       double estimate = fw::estimate_field(*(iData.innerTrack()), true);
0217       if (estimate >= 0)
0218         field->guessField(estimate);
0219     }
0220     if (iData.outerTrack().isAvailable()) {
0221       double estimate = fw::estimate_field(*(iData.outerTrack()));
0222       if (estimate >= 0)
0223         field->guessFieldIsOn(estimate > 0.5);
0224     }
0225   }
0226 }
0227
0228 //______________________________________________________________________________
0229
0230 void FWMuonBuilder::buildMuon(
0231     FWProxyBuilderBase* pb, const reco::Muon* muon, TEveElement* tList, bool showEndcap, bool tracksOnly) {
0232   calculateField(*muon, pb->context().getField());
0233
0234   TEveRecTrack recTrack;
0235   recTrack.fBeta = 1.;
0236
0237   // If we deal with a tracker muon we use the inner track and guide it
0238   // through the trajectory points from the reconstruction. Segments
0239   // represent hits. Matching between hits and the trajectory shows
0240   // how well the inner track matches with the muon hypothesis.
0241   //
0242   // In other cases we use a global muon track with a few states from
0243   // the inner and outer tracks or just the outer track if it's the
0244   // only option
0245
0246   if (muon->isTrackerMuon() && muon->innerTrack().isAvailable() && muon->isMatchesValid() &&
0247       !buggyMuon(&*muon, pb->context().getGeom())) {
0248     TEveTrack* trk = fireworks::prepareTrack(
0249         *(muon->innerTrack()), pb->context().getMuonTrackPropagator(), getRecoTrajectoryPoints(muon, pb->item()));
0250     trk->MakeTrack();
0251     trk->SetLineWidth(m_lineWidth);
0252     pb->setupAddElement(trk, tList);
0253     if (!tracksOnly)
0254       addMatchInformation(&(*muon), pb, tList, showEndcap);
0255     return;
0256   }
0257
0258   if (muon->isGlobalMuon() && muon->globalTrack().isAvailable()) {
0259     std::vector<TEveVector> extraPoints;
0260     if (muon->innerTrack().isAvailable() && muon->innerTrack()->extra().isAvailable()) {
0261       extraPoints.push_back(TEveVector(muon->innerTrack()->innerPosition().x(),
0262                                        muon->innerTrack()->innerPosition().y(),
0263                                        muon->innerTrack()->innerPosition().z()));
0264       extraPoints.push_back(TEveVector(muon->innerTrack()->outerPosition().x(),
0265                                        muon->innerTrack()->outerPosition().y(),
0266                                        muon->innerTrack()->outerPosition().z()));
0267     }
0268     if (muon->outerTrack().isAvailable() && muon->outerTrack()->extra().isAvailable()) {
0269       extraPoints.push_back(TEveVector(muon->outerTrack()->innerPosition().x(),
0270                                        muon->outerTrack()->innerPosition().y(),
0271                                        muon->outerTrack()->innerPosition().z()));
0272       extraPoints.push_back(TEveVector(muon->outerTrack()->outerPosition().x(),
0273                                        muon->outerTrack()->outerPosition().y(),
0274                                        muon->outerTrack()->outerPosition().z()));
0275     }
0276     TEveTrack* trk = nullptr;
0277     if (extraPoints.empty())
0278       trk = fireworks::prepareTrack(*(muon->globalTrack()), pb->context().getMuonTrackPropagator());
0279     else
0280       trk =
0281           prepareMuonTrackWithExtraPoints(*(muon->globalTrack()), pb->context().getMuonTrackPropagator(), extraPoints);
0282
0283     trk->MakeTrack();
0284     trk->SetLineWidth(m_lineWidth);
0285     pb->setupAddElement(trk, tList);
0286     return;
0287   }
0288
0289   if (muon->innerTrack().isAvailable()) {
0290     TEveTrack* trk = fireworks::prepareTrack(*(muon->innerTrack()), pb->context().getMuonTrackPropagator());
0291     trk->MakeTrack();
0292     pb->setupAddElement(trk, tList);
0293     return;
0294   }
0295
0296   if (muon->outerTrack().isAvailable()) {
0297     TEveTrack* trk = fireworks::prepareTrack(*(muon->outerTrack()), pb->context().getMuonTrackPropagator());
0298     trk->MakeTrack();
0299     trk->SetLineWidth(m_lineWidth);
0300     pb->setupAddElement(trk, tList);
0301     return;
0302   }
0303
0304   // if got that far it means we have nothing but a candidate
0305   // show it anyway.
0306   TEveTrack* trk = fireworks::prepareCandidate(*muon, pb->context().getMuonTrackPropagator());
0307   trk->MakeTrack();
0308   trk->SetLineWidth(m_lineWidth);
0309   pb->setupAddElement(trk, tList);
0310 }