Tracks/src/TrackUtils.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     Tracks
0004 // Class  :     TrackUtils
0005 //
0006
0007 // system include files
0008 #include "TEveTrack.h"
0009 #include "TEveTrackPropagator.h"
0010 #include "TEveStraightLineSet.h"
0011 #include "TEveVSDStructs.h"
0012 #include "TEveGeoNode.h"
0013
0014 // user include files
0015 #include "Fireworks/Tracks/interface/TrackUtils.h"
0016 #include "Fireworks/Core/interface/FWEventItem.h"
0017 #include "Fireworks/Core/interface/FWGeometry.h"
0018 #include "Fireworks/Core/interface/TEveElementIter.h"
0019 #include "Fireworks/Core/interface/fwLog.h"
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024
0025 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0026
0027 #include "DataFormats/TrackingRecHit/interface/RecHit2DLocalPos.h"
0028 #include "DataFormats/TrackingRecHit/interface/RecSegment.h"
0029 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0030 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0031 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0033 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0034 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0035 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0036
0037 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0038 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0039 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0040 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0041 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0042 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0043
0044 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0045 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0046 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0047
0048 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0049
0050 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0051 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0052
0053 #include "FWCore/Common/interface/EventBase.h"
0054
0055 namespace fireworks {
0056
0057   static const double MICRON = 1. / 1000. / 1000.;
0058
0059   // -- Si module names for printout
0060   static const std::string subdets[7] = {"UNKNOWN", "PXB", "PXF", "TIB", "TID", "TOB", "TEC"};
0061
0062   TEveTrack* prepareTrack(const reco::Track& track,
0063                           TEveTrackPropagator* propagator,
0064                           const std::vector<TEveVector>& extraRefPoints) {
0065     // To make use of all available information, we have to order states
0066     // properly first. Propagator should take care of y=0 transition.
0067
0068     std::vector<State> refStates;
0069     TEveVector trackMomentum(track.px(), track.py(), track.pz());
0070     refStates.push_back(State(TEveVector(track.vx(), track.vy(), track.vz()), trackMomentum));
0071     if (track.extra().isAvailable()) {
0072       if (track.innerOk()) {
0073         const reco::TrackBase::Point& v = track.innerPosition();
0074         const reco::TrackBase::Vector& p = track.innerMomentum();
0075         refStates.push_back(State(TEveVector(v.x(), v.y(), v.z()), TEveVector(p.x(), p.y(), p.z())));
0076       }
0077       if (track.outerOk()) {
0078         const reco::TrackBase::Point& v = track.outerPosition();
0079         const reco::TrackBase::Vector& p = track.outerMomentum();
0080         refStates.push_back(State(TEveVector(v.x(), v.y(), v.z()), TEveVector(p.x(), p.y(), p.z())));
0081       }
0082     }
0083     for (std::vector<TEveVector>::const_iterator point = extraRefPoints.begin(), pointEnd = extraRefPoints.end();
0084          point != pointEnd;
0085          ++point)
0086       refStates.push_back(State(*point));
0087     if (track.pt() > 1)
0088       std::sort(refStates.begin(), refStates.end(), StateOrdering(trackMomentum));
0089
0090     // * if the first state has non-zero momentum use it as a starting point
0091     //   and all other points as PathMarks to follow
0092     // * if the first state has only position, try the last state. If it has
0093     //   momentum we propagate backword, if not, we look for the first one
0094     //   on left that has momentum and ignore all earlier.
0095     //
0096
0097     TEveRecTrack t;
0098     t.fBeta = 1.;
0099     t.fSign = track.charge();
0100
0101     if (refStates.front().valid) {
0102       t.fV = refStates.front().position;
0103       t.fP = refStates.front().momentum;
0104       TEveTrack* trk = new TEveTrack(&t, propagator);
0105       for (unsigned int i(1); i < refStates.size() - 1; ++i) {
0106         if (refStates[i].valid)
0107           trk->AddPathMark(TEvePathMark(TEvePathMark::kReference, refStates[i].position, refStates[i].momentum));
0108         else
0109           trk->AddPathMark(TEvePathMark(TEvePathMark::kDaughter, refStates[i].position));
0110       }
0111       if (refStates.size() > 1) {
0112         trk->AddPathMark(TEvePathMark(TEvePathMark::kDecay, refStates.back().position));
0113       }
0114       return trk;
0115     }
0116
0117     if (refStates.back().valid) {
0118       t.fSign = (-1) * track.charge();
0119       t.fV = refStates.back().position;
0120       t.fP = refStates.back().momentum * (-1.0f);
0121       TEveTrack* trk = new TEveTrack(&t, propagator);
0122       unsigned int i(refStates.size() - 1);
0123       for (; i > 0; --i) {
0124         if (refStates[i].valid)
0125           trk->AddPathMark(
0126               TEvePathMark(TEvePathMark::kReference, refStates[i].position, refStates[i].momentum * (-1.0f)));
0127         else
0128           trk->AddPathMark(TEvePathMark(TEvePathMark::kDaughter, refStates[i].position));
0129       }
0130       if (refStates.size() > 1) {
0131         trk->AddPathMark(TEvePathMark(TEvePathMark::kDecay, refStates.front().position));
0132       }
0133       return trk;
0134     }
0135
0136     unsigned int i(0);
0137     while (i < refStates.size() && !refStates[i].valid)
0138       ++i;
0139     assert(i < refStates.size());
0140
0141     t.fV = refStates[i].position;
0142     t.fP = refStates[i].momentum;
0143     TEveTrack* trk = new TEveTrack(&t, propagator);
0144     for (unsigned int j(i + 1); j < refStates.size() - 1; ++j) {
0145       if (refStates[i].valid)
0146         trk->AddPathMark(TEvePathMark(TEvePathMark::kReference, refStates[i].position, refStates[i].momentum));
0147       else
0148         trk->AddPathMark(TEvePathMark(TEvePathMark::kDaughter, refStates[i].position));
0149     }
0150     if (i < refStates.size()) {
0151       trk->AddPathMark(TEvePathMark(TEvePathMark::kDecay, refStates.back().position));
0152     }
0153     return trk;
0154   }
0155
0156   //==============================================================================
0157
0158   // Transform measurement to local coordinates in X dimension
0159   float pixelLocalX(const double mpx, const float* par) {
0160     float xoffset = 0;
0161     float xpitch = 0;
0162     // std::cerr << "version pr0d " << fireworks::Context::getInstance()->getGeom()->getProducerVersion() << "\n";
0163     if (fireworks::Context::getInstance() && fireworks::Context::getInstance()->getGeom()->getProducerVersion() < 1) {
0164       static const int ROWS_PER_ROC = 80;      // Num of cols per ROC
0165       static const int BIG_PIX_PER_ROC_X = 1;  // in x direction, rows
0166       static const double PITCHX = 100 * MICRON;
0167       // Calculate the edge of the active sensor with respect to the center,
0168       // that is simply the half-size.
0169       // Take into account large pixels
0170       xpitch = PITCHX;
0171       xoffset = -(par[0] + BIG_PIX_PER_ROC_X * par[0] / ROWS_PER_ROC) / 2. * PITCHX;
0172     } else {
0173       xpitch = par[0];
0174       xoffset = par[2];
0175     }
0176
0177     bool bigPixelsLayout = par[4];
0178
0179     int binoffx = int(mpx);            // truncate to int
0180     double fractionX = mpx - binoffx;  // find the fraction
0181     double local_xpitch = xpitch;      // defaultpitch
0182
0183     if (bigPixelsLayout == 0) {
0184       // Measurement to local transformation for X coordinate
0185       // X coordinate is in the ROC row number direction
0186       // Copy from RectangularPixelTopology::localX implementation
0187       if (binoffx > 80) {  // ROC 1 - handles x on edge cluster
0188         binoffx = binoffx + 2;
0189       } else if (binoffx == 80) {  // ROC 1
0190         binoffx = binoffx + 1;
0191         local_xpitch = 2 * xpitch;
0192       } else if (binoffx == 79) {  // ROC 0
0193         binoffx = binoffx + 0;
0194         local_xpitch = 2 * xpitch;
0195       } else if (binoffx >= 0) {  // ROC 0
0196         binoffx = binoffx + 0;
0197       }
0198     }
0199
0200     // The final position in local coordinates
0201     double lpX = double(binoffx * xpitch) + fractionX * local_xpitch + xoffset;
0202
0203     return lpX;
0204   }
0205
0206   //==============================================================================
0207
0208   // Transform measurement to local coordinates in Y dimension
0209   float pixelLocalY(const double mpy, const float* par) {
0210     float ypitch = 0;
0211     float yoffset = 0;
0212     if (fireworks::Context::getInstance() && fireworks::Context::getInstance()->getGeom()->getProducerVersion() < 1) {
0213       static const double PITCHY = 150 * MICRON;
0214       static const int BIG_PIX_PER_ROC_Y = 2;  // in y direction, cols
0215       static const int COLS_PER_ROC = 52;      // Num of Rows per ROC
0216
0217       // Calculate the edge of the active sensor with respect to the center,
0218       // that is simply the half-size.
0219       // Take into account large pixels
0220       yoffset = -(par[1] + BIG_PIX_PER_ROC_Y * par[1] / COLS_PER_ROC) / 2. * PITCHY;
0221       ypitch = PITCHY;
0222     } else {
0223       ypitch = par[1];
0224       yoffset = par[3];
0225     }
0226     // Measurement to local transformation for Y coordinate
0227     // Y is in the ROC column number direction
0228     // Copy from RectangularPixelTopology::localY implementation
0229     int binoffy = int(mpy);            // truncate to int
0230     double fractionY = mpy - binoffy;  // find the fraction
0231     double local_pitchy = ypitch;      // defaultpitch
0232
0233     bool bigPixelsLayout = par[4];
0234     if (bigPixelsLayout == 0) {
0235       constexpr int bigYIndeces[]{0, 51, 52, 103, 104, 155, 156, 207, 208, 259, 260, 311, 312, 363, 364, 415, 416, 511};
0236       auto const j = std::lower_bound(std::begin(bigYIndeces), std::end(bigYIndeces), binoffy);
0237       if (*j == binoffy)
0238         local_pitchy *= 2;
0239       binoffy += (j - bigYIndeces);
0240     }
0241
0242     // The final position in local coordinates
0243     double lpY = double(binoffy * ypitch) + fractionY * local_pitchy + yoffset;
0244
0245     return lpY;
0246   }
0247
0248   // Transform measurement to local coordinates in X dimension
0249   float phase2PixelLocalX(const double mpx, const float* par, const float* shape) {
0250     // The final position in local coordinates
0251     return (-shape[1] + mpx * par[0]);
0252   }
0253
0254   // Transform measurement to local coordinates in Y dimension
0255   float phase2PixelLocalY(const double mpy, const float* par, const float* shape) {
0256     // The final position in local coordinates
0257     return (-shape[2] + mpy * par[1]);
0258   }
0259
0260   //
0261   // Returns strip geometry in local coordinates of a detunit.
0262   // The strip is a line from a localTop to a localBottom point.
0263   void localSiStrip(short strip, float* localTop, float* localBottom, const float* pars, unsigned int id) {
0264     Float_t topology = pars[0];
0265     Float_t halfStripLength = pars[2] * 0.5;
0266
0267     Double_t localCenter[3] = {0.0, 0.0, 0.0};
0268     localTop[1] = halfStripLength;
0269     localBottom[1] = -halfStripLength;
0270
0271     if (topology == 1)  // RadialStripTopology
0272     {
0273       // stripAngle = phiOfOneEdge + strip * angularWidth
0274       // localY = originToIntersection * tan( stripAngle )
0275       Float_t stripAngle = tan(pars[5] + strip * pars[6]);
0276       Float_t delta = halfStripLength * stripAngle;
0277       localCenter[0] = pars[4] * stripAngle;
0278       localTop[0] = localCenter[0] + delta;
0279       localBottom[0] = localCenter[0] - delta;
0280     } else if (topology == 2)  // RectangularStripTopology
0281     {
0282       // offset = -numberOfStrips/2. * pitch
0283       // localY = strip * pitch + offset
0284       Float_t offset = -pars[1] * 0.5 * pars[3];
0285       localCenter[0] = strip * pars[3] + offset;
0286       localTop[0] = localCenter[0];
0287       localBottom[0] = localCenter[0];
0288     } else if (topology == 3)  // TrapezoidalStripTopology
0289     {
0290       fwLog(fwlog::kError) << "did not expect TrapezoidalStripTopology of " << id << std::endl;
0291     } else if (pars[0] == 0)  // StripTopology
0292     {
0293       fwLog(fwlog::kError) << "did not find StripTopology of " << id << std::endl;
0294     }
0295   }
0296
0297   //______________________________________________________________________________
0298
0299   void setupAddElement(TEveElement* el, TEveElement* parent, const FWEventItem* item, bool master, bool color) {
0300     if (master) {
0301       el->CSCTakeAnyParentAsMaster();
0302       el->SetPickable(true);
0303     }
0304
0305     if (color) {
0306       el->CSCApplyMainColorToMatchingChildren();
0307       el->CSCApplyMainTransparencyToMatchingChildren();
0308       el->SetMainColor(item->defaultDisplayProperties().color());
0309       assert((item->defaultDisplayProperties().transparency() >= 0) &&
0310              (item->defaultDisplayProperties().transparency() <= 100));
0311       el->SetMainTransparency(item->defaultDisplayProperties().transparency());
0312     }
0313     parent->AddElement(el);
0314   }
0315
0316   //______________________________________________________________________________
0317
0318   const SiStripCluster* extractClusterFromTrackingRecHit(const TrackingRecHit* rechit) {
0319     const SiStripCluster* cluster = nullptr;
0320
0321     if (const SiStripRecHit2D* hit2D = dynamic_cast<const SiStripRecHit2D*>(rechit)) {
0322       fwLog(fwlog::kDebug) << "hit 2D ";
0323
0324       cluster = hit2D->cluster().get();
0325     }
0326     if (cluster == nullptr) {
0327       if (const SiStripRecHit1D* hit1D = dynamic_cast<const SiStripRecHit1D*>(rechit)) {
0328         fwLog(fwlog::kDebug) << "hit 1D ";
0329
0330         cluster = hit1D->cluster().get();
0331       }
0332     }
0333     return cluster;
0334   }
0335
0336   void addSiStripClusters(
0337       const FWEventItem* iItem, const reco::Track& t, class TEveElement* tList, bool addNearbyClusters, bool master) {
0338     // master is true if the product is for proxy builder
0339     const FWGeometry* geom = iItem->getGeom();
0340
0341     const edmNew::DetSetVector<SiStripCluster>* allClusters = nullptr;
0342     if (addNearbyClusters) {
0343       for (trackingRecHit_iterator it = t.recHitsBegin(), itEnd = t.recHitsEnd(); it != itEnd; ++it) {
0344         const auto& rhs = *(*(it));
0345         if (typeid(rhs) == typeid(SiStripRecHit2D)) {
0346           const SiStripRecHit2D& hit = static_cast<const SiStripRecHit2D&>(**it);
0347           if (hit.cluster().isNonnull() && hit.cluster().isAvailable()) {
0348             edm::Handle<edmNew::DetSetVector<SiStripCluster> > allClustersHandle;
0349             iItem->getEvent()->get(hit.cluster().id(), allClustersHandle);
0350             allClusters = allClustersHandle.product();
0351             break;
0352           }
0353         } else if (typeid(rhs) == typeid(SiStripRecHit1D)) {
0354           const SiStripRecHit1D& hit = static_cast<const SiStripRecHit1D&>(**it);
0355           if (hit.cluster().isNonnull() && hit.cluster().isAvailable()) {
0356             edm::Handle<edmNew::DetSetVector<SiStripCluster> > allClustersHandle;
0357             iItem->getEvent()->get(hit.cluster().id(), allClustersHandle);
0358             allClusters = allClustersHandle.product();
0359             break;
0360           }
0361         }
0362       }
0363     }
0364
0365     for (trackingRecHit_iterator it = t.recHitsBegin(), itEnd = t.recHitsEnd(); it != itEnd; ++it) {
0366       unsigned int rawid = (*it)->geographicalId();
0367       if (!geom->contains(rawid)) {
0368         fwLog(fwlog::kError) << "failed to get geometry of SiStripCluster with detid: " << rawid << std::endl;
0369
0370         continue;
0371       }
0372
0373       const float* pars = geom->getParameters(rawid);
0374
0375       // -- get phi from SiStripHit
0376       auto rechitRef = *it;
0377       const TrackingRecHit* rechit = &(*rechitRef);
0378       const SiStripCluster* cluster = extractClusterFromTrackingRecHit(rechit);
0379
0380       if (cluster) {
0381         if (allClusters != nullptr) {
0382           const edmNew::DetSet<SiStripCluster>& clustersOnThisDet = (*allClusters)[rechit->geographicalId().rawId()];
0383
0384           for (edmNew::DetSet<SiStripCluster>::const_iterator itc = clustersOnThisDet.begin(),
0385                                                               edc = clustersOnThisDet.end();
0386                itc != edc;
0387                ++itc) {
0388             TEveStraightLineSet* scposition = new TEveStraightLineSet;
0389             scposition->SetDepthTest(false);
0390             scposition->SetPickable(kTRUE);
0391
0392             short firststrip = itc->firstStrip();
0393
0394             if (&*itc == cluster) {
0395               scposition->SetTitle(Form("Exact SiStripCluster from TrackingRecHit, first strip %d", firststrip));
0396               scposition->SetLineColor(kGreen);
0397             } else {
0398               scposition->SetTitle(Form("SiStripCluster, first strip %d", firststrip));
0399               scposition->SetLineColor(kRed);
0400             }
0401
0402             float localTop[3] = {0.0, 0.0, 0.0};
0403             float localBottom[3] = {0.0, 0.0, 0.0};
0404
0405             fireworks::localSiStrip(firststrip, localTop, localBottom, pars, rawid);
0406
0407             float globalTop[3];
0408             float globalBottom[3];
0409             geom->localToGlobal(rawid, localTop, globalTop, localBottom, globalBottom);
0410
0411             scposition->AddLine(
0412                 globalTop[0], globalTop[1], globalTop[2], globalBottom[0], globalBottom[1], globalBottom[2]);
0413
0414             setupAddElement(scposition, tList, iItem, master, false);
0415           }
0416         } else {
0417           short firststrip = cluster->firstStrip();
0418           TEveStraightLineSet* scposition = new TEveStraightLineSet;
0419           scposition->SetDepthTest(false);
0420           scposition->SetPickable(kTRUE);
0421           scposition->SetTitle(Form("SiStripCluster, first strip %d", firststrip));
0422
0423           float localTop[3] = {0.0, 0.0, 0.0};
0424           float localBottom[3] = {0.0, 0.0, 0.0};
0425
0426           fireworks::localSiStrip(firststrip, localTop, localBottom, pars, rawid);
0427
0428           float globalTop[3];
0429           float globalBottom[3];
0430           geom->localToGlobal(rawid, localTop, globalTop, localBottom, globalBottom);
0431
0432           scposition->AddLine(
0433               globalTop[0], globalTop[1], globalTop[2], globalBottom[0], globalBottom[1], globalBottom[2]);
0434
0435           setupAddElement(scposition, tList, iItem, master, true);
0436         }
0437       } else if (!rechit->isValid() && (rawid != 0))  // lost hit
0438       {
0439         if (allClusters != nullptr) {
0440           edmNew::DetSetVector<SiStripCluster>::const_iterator itds = allClusters->find(rawid);
0441           if (itds != allClusters->end()) {
0442             const edmNew::DetSet<SiStripCluster>& clustersOnThisDet = *itds;
0443             for (edmNew::DetSet<SiStripCluster>::const_iterator itc = clustersOnThisDet.begin(),
0444                                                                 edc = clustersOnThisDet.end();
0445                  itc != edc;
0446                  ++itc) {
0447               short firststrip = itc->firstStrip();
0448
0449               TEveStraightLineSet* scposition = new TEveStraightLineSet;
0450               scposition->SetDepthTest(false);
0451               scposition->SetPickable(kTRUE);
0452               scposition->SetTitle(Form("Lost SiStripCluster, first strip %d", firststrip));
0453
0454               float localTop[3] = {0.0, 0.0, 0.0};
0455               float localBottom[3] = {0.0, 0.0, 0.0};
0456
0457               fireworks::localSiStrip(firststrip, localTop, localBottom, pars, rawid);
0458
0459               float globalTop[3];
0460               float globalBottom[3];
0461               geom->localToGlobal(rawid, localTop, globalTop, localBottom, globalBottom);
0462
0463               scposition->AddLine(
0464                   globalTop[0], globalTop[1], globalTop[2], globalBottom[0], globalBottom[1], globalBottom[2]);
0465
0466               setupAddElement(scposition, tList, iItem, master, false);
0467               scposition->SetLineColor(kRed);
0468             }
0469           }
0470         }
0471       } else {
0472         fwLog(fwlog::kDebug) << "*ANOTHER* option possible: valid=" << rechit->isValid() << ", rawid=" << rawid
0473                              << std::endl;
0474       }
0475     }
0476   }
0477
0478   //______________________________________________________________________________
0479
0480   void pushNearbyPixelHits(std::vector<TVector3>& pixelPoints, const FWEventItem& iItem, const reco::Track& t) {
0481     const edmNew::DetSetVector<SiPixelCluster>* allClusters = nullptr;
0482     for (trackingRecHit_iterator it = t.recHitsBegin(), itEnd = t.recHitsEnd(); it != itEnd; ++it) {
0483       const auto& rhs = *(*(it));
0484       if (typeid(rhs) == typeid(SiPixelRecHit)) {
0485         const SiPixelRecHit& hit = static_cast<const SiPixelRecHit&>(**it);
0486         if (hit.cluster().isNonnull() && hit.cluster().isAvailable()) {
0487           edm::Handle<edmNew::DetSetVector<SiPixelCluster> > allClustersHandle;
0488           iItem.getEvent()->get(hit.cluster().id(), allClustersHandle);
0489           allClusters = allClustersHandle.product();
0490           break;
0491         }
0492       }
0493     }
0494     if (allClusters == nullptr)
0495       return;
0496
0497     const FWGeometry* geom = iItem.getGeom();
0498
0499     for (trackingRecHit_iterator it = t.recHitsBegin(), itEnd = t.recHitsEnd(); it != itEnd; ++it) {
0500       const TrackingRecHit* rh = &(**it);
0501
0502       DetId id = (*it)->geographicalId();
0503       if (!geom->contains(id)) {
0504         fwLog(fwlog::kError) << "failed to get geometry of Tracker Det with raw id: " << id.rawId() << std::endl;
0505
0506         continue;
0507       }
0508
0509       // -- in which detector are we?
0510       unsigned int subdet = (unsigned int)id.subdetId();
0511       if ((subdet != PixelSubdetector::PixelBarrel) && (subdet != PixelSubdetector::PixelEndcap))
0512         continue;
0513
0514       const SiPixelCluster* hitCluster = nullptr;
0515       if (const SiPixelRecHit* pixel = dynamic_cast<const SiPixelRecHit*>(rh))
0516         hitCluster = pixel->cluster().get();
0517       edmNew::DetSetVector<SiPixelCluster>::const_iterator itds = allClusters->find(id.rawId());
0518       if (itds != allClusters->end()) {
0519         const edmNew::DetSet<SiPixelCluster>& clustersOnThisDet = *itds;
0520         for (edmNew::DetSet<SiPixelCluster>::const_iterator itc = clustersOnThisDet.begin(),
0521                                                             edc = clustersOnThisDet.end();
0522              itc != edc;
0523              ++itc) {
0524           if (&*itc != hitCluster)
0525             pushPixelCluster(pixelPoints, *geom, id, *itc, geom->getParameters(id));
0526         }
0527       }
0528     }
0529   }
0530
0531   //______________________________________________________________________________
0532
0533   void pushPixelHits(std::vector<TVector3>& pixelPoints, const FWEventItem& iItem, const reco::Track& t) {
0534     /*
0535     * -- return for each Pixel Hit a 3D point
0536     */
0537     const FWGeometry* geom = iItem.getGeom();
0538
0539     double dz = t.dz();
0540     double vz = t.vz();
0541     double etaT = t.eta();
0542
0543     fwLog(fwlog::kDebug) << "Track eta: " << etaT << ", vz: " << vz << ", dz: " << dz << std::endl;
0544
0545     int cnt = 0;
0546     for (trackingRecHit_iterator it = t.recHitsBegin(), itEnd = t.recHitsEnd(); it != itEnd; ++it) {
0547       const TrackingRecHit* rh = &(**it);
0548       // -- get position of center of wafer, assuming (0,0,0) is the center
0549       DetId id = (*it)->geographicalId();
0550       if (!geom->contains(id)) {
0551         fwLog(fwlog::kError) << "failed to get geometry of Tracker Det with raw id: " << id.rawId() << std::endl;
0552
0553         continue;
0554       }
0555
0556       // -- in which detector are we?
0557       unsigned int subdet = (unsigned int)id.subdetId();
0558
0559       if ((subdet == PixelSubdetector::PixelBarrel) || (subdet == PixelSubdetector::PixelEndcap)) {
0560         fwLog(fwlog::kDebug) << cnt++ << " -- " << subdets[subdet];
0561
0562         if (const SiPixelRecHit* pixel = dynamic_cast<const SiPixelRecHit*>(rh)) {
0563           const SiPixelCluster& c = *(pixel->cluster());
0564           pushPixelCluster(pixelPoints, *geom, id, c, geom->getParameters(id));
0565         }
0566       }
0567     }
0568   }
0569
0570   void pushPixelCluster(std::vector<TVector3>& pixelPoints,
0571                         const FWGeometry& geom,
0572                         DetId id,
0573                         const SiPixelCluster& c,
0574                         const float* pars) {
0575     double row = c.minPixelRow();
0576     double col = c.minPixelCol();
0577     float lx = 0.;
0578     float ly = 0.;
0579
0580     // int nrows = (int)pars[0];
0581     // int ncols = (int)pars[1];
0582     lx = pixelLocalX(row, pars);
0583     ly = pixelLocalY(col, pars);
0584
0585     fwLog(fwlog::kDebug) << ", row: " << row << ", col: " << col << ", lx: " << lx << ", ly: " << ly;
0586
0587     float local[3] = {lx, ly, 0.};
0588     float global[3];
0589     geom.localToGlobal(id, local, global);
0590     TVector3 pb(global[0], global[1], global[2]);
0591     pixelPoints.push_back(pb);
0592
0593     fwLog(fwlog::kDebug) << " x: " << pb.X() << ", y: " << pb.Y() << " z: " << pb.Z() << " eta: " << pb.Eta()
0594                          << ", phi: " << pb.Phi() << " rho: " << pb.Pt() << std::endl;
0595   }
0596
0597   //______________________________________________________________________________
0598
0599   std::string info(const DetId& id) {
0600     std::ostringstream oss;
0601
0602     oss << "DetId: " << id.rawId() << "\n";
0603
0604     switch (id.det()) {
0605       case DetId::Tracker:
0606         oss << fireworks::Context::getInstance()->getGeom()->getTrackerTopology()->print(id.rawId());
0607         break;
0608
0609       case DetId::Muon:
0610         switch (id.subdetId()) {
0611           case MuonSubdetId::DT: {
0612             DTChamberId detId(id.rawId());
0613             oss << "DT chamber (wheel, station, sector): " << detId.wheel() << ", " << detId.station() << ", "
0614                 << detId.sector();
0615           } break;
0616           case MuonSubdetId::CSC: {
0617             CSCDetId detId(id.rawId());
0618             oss << "CSC chamber (endcap, station, ring, chamber, layer): " << detId.endcap() << ", " << detId.station()
0619                 << ", " << detId.ring() << ", " << detId.chamber() << ", " << detId.layer();
0620           } break;
0621           case MuonSubdetId::RPC: {
0622             RPCDetId detId(id.rawId());
0623             oss << "RPC chamber ";
0624             switch (detId.region()) {
0625               case 0:
0626                 oss << "/ barrel / (wheel, station, sector, layer, subsector, roll): " << detId.ring() << ", "
0627                     << detId.station() << ", " << detId.sector() << ", " << detId.layer() << ", " << detId.subsector()
0628                     << ", " << detId.roll();
0629                 break;
0630               case 1:
0631                 oss << "/ forward endcap / (wheel, station, sector, layer, subsector, roll): " << detId.ring() << ", "
0632                     << detId.station() << ", " << detId.sector() << ", " << detId.layer() << ", " << detId.subsector()
0633                     << ", " << detId.roll();
0634                 break;
0635               case -1:
0636                 oss << "/ backward endcap / (wheel, station, sector, layer, subsector, roll): " << detId.ring() << ", "
0637                     << detId.station() << ", " << detId.sector() << ", " << detId.layer() << ", " << detId.subsector()
0638                     << ", " << detId.roll();
0639                 break;
0640             }
0641           } break;
0642           case MuonSubdetId::GEM: {
0643             GEMDetId detId(id.rawId());
0644             oss << "GEM chamber (region, station, ring, chamber, layer): " << detId.region() << ", " << detId.station()
0645                 << ", " << detId.ring() << ", " << detId.chamber() << ", " << detId.layer();
0646           } break;
0647           case MuonSubdetId::ME0: {
0648             ME0DetId detId(id.rawId());
0649             oss << "ME0 chamber (region, chamber, layer): " << detId.region() << ", " << detId.chamber() << ", "
0650                 << detId.layer();
0651           } break;
0652         }
0653         break;
0654
0655       case DetId::Calo: {
0656         CaloTowerDetId detId(id.rawId());
0657         oss << "CaloTower (ieta, iphi): " << detId.ieta() << ", " << detId.iphi();
0658       } break;
0659
0660       case DetId::Ecal:
0661         switch (id.subdetId()) {
0662           case EcalBarrel: {
0663             EBDetId detId(id);
0664             oss << "EcalBarrel (ieta, iphi, tower_ieta, tower_iphi): " << detId.ieta() << ", " << detId.iphi() << ", "
0665                 << detId.tower_ieta() << ", " << detId.tower_iphi();
0666           } break;
0667           case EcalEndcap: {
0668             EEDetId detId(id);
0669             oss << "EcalEndcap (ix, iy, SuperCrystal, crystal, quadrant): " << detId.ix() << ", " << detId.iy() << ", "
0670                 << detId.isc() << ", " << detId.ic() << ", " << detId.iquadrant();
0671           } break;
0672           case EcalPreshower:
0673             oss << "EcalPreshower";
0674             break;
0675           case EcalTriggerTower:
0676             oss << "EcalTriggerTower";
0677             break;
0678           case EcalLaserPnDiode:
0679             oss << "EcalLaserPnDiode";
0680             break;
0681         }
0682         break;
0683
0684       case DetId::Hcal: {
0685         HcalDetId detId(id);
0686         switch (detId.subdet()) {
0687           case HcalEmpty:
0688             oss << "HcalEmpty ";
0689             break;
0690           case HcalBarrel:
0691             oss << "HcalBarrel ";
0692             break;
0693           case HcalEndcap:
0694             oss << "HcalEndcap ";
0695             break;
0696           case HcalOuter:
0697             oss << "HcalOuter ";
0698             break;
0699           case HcalForward:
0700             oss << "HcalForward ";
0701             break;
0702           case HcalTriggerTower:
0703             oss << "HcalTriggerTower ";
0704             break;
0705           case HcalOther:
0706             oss << "HcalOther ";
0707             break;
0708         }
0709         oss << "(ieta, iphi, depth):" << detId.ieta() << ", " << detId.iphi() << ", " << detId.depth();
0710       } break;
0711       default:;
0712     }
0713     return oss.str();
0714   }
0715
0716   std::string info(const std::set<DetId>& idSet) {
0717     std::string text;
0718     for (std::set<DetId>::const_iterator id = idSet.begin(), idEnd = idSet.end(); id != idEnd; ++id) {
0719       text += info(*id);
0720       text += "\n";
0721     }
0722     return text;
0723   }
0724
0725   std::string info(const std::vector<DetId>& idSet) {
0726     std::string text;
0727     for (std::vector<DetId>::const_iterator id = idSet.begin(), idEnd = idSet.end(); id != idEnd; ++id) {
0728       text += info(*id);
0729       text += "\n";
0730     }
0731     return text;
0732   }
0733 }  // namespace fireworks