ParticleFlow/plugins/FWPFCandidateWithHitsProxyBuilder.cc

0001 #include "TEveBoxSet.h"
0002 #include "TEveTrack.h"
0003 #include "TEveTrackPropagator.h"
0004 #include "TEveCompound.h"
0005 #include "TEveStraightLineSet.h"
0006 #include "TEveProjectionBases.h"
0007
0008 #include "Fireworks/ParticleFlow/plugins/FWPFCandidateWithHitsProxyBuilder.h"
0009 #include "Fireworks/Core/interface/FWEventItem.h"
0010 #include "Fireworks/Core/interface/FWGeometry.h"
0011 #include "Fireworks/Core/interface/BuilderUtils.h"
0012 #include "Fireworks/Core/interface/fwLog.h"
0013 #include "Fireworks/Core/interface/FWViewEnergyScale.h"
0014 #include "Fireworks/ParticleFlow/interface/setTrackTypePF.h"
0015
0016 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0017 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0020
0021 #include "DataFormats/FWLite/interface/Handle.h"
0022
0023 namespace {
0024   const std::string cname("particleFlowRecHitHCALUpgrade");
0025
0026   void addLineToLineSet(TEveStraightLineSet* ls, const float* p, int i1, int i2) {
0027     i1 *= 3;
0028     i2 *= 3;
0029     ls->AddLine(p[i1], p[i1 + 1], p[i1 + 2], p[i2], p[i2 + 1], p[i2 + 2]);
0030   }
0031
0032   void addBoxAsLines(TEveStraightLineSet* lineset, const float* p) {
0033     for (int l = 0; l < 5; l += 4) {
0034       addLineToLineSet(lineset, p, 0 + l, 1 + l);
0035       addLineToLineSet(lineset, p, 1 + l, 2 + l);
0036       addLineToLineSet(lineset, p, 2 + l, 3 + l);
0037       addLineToLineSet(lineset, p, 3 + l, 0 + l);
0038     }
0039     for (int l = 0; l < 4; ++l)
0040       addLineToLineSet(lineset, p, 0 + l, 4 + l);
0041   }
0042
0043   void editLineInLineSet(TEveChunkManager::iterator& li, const float* p, int i1, int i2) {
0044     TEveStraightLineSet::Line_t& line = *(TEveStraightLineSet::Line_t*)li();
0045     i1 *= 3;
0046     i2 *= 3;
0047     for (int i = 0; i < 3; ++i) {
0048       line.fV1[0 + i] = p[i1 + i];
0049       line.fV2[0 + i] = p[i2 + i];
0050     }
0051
0052     li.next();
0053   }
0054
0055   void editBoxInLineSet(TEveChunkManager::iterator& li, const float* p) {
0056     for (int i = 0; i < 5; i += 4) {
0057       editLineInLineSet(li, p, 0 + i, 1 + i);
0058
0059       editLineInLineSet(li, p, 1 + i, 2 + i);
0060       editLineInLineSet(li, p, 2 + i, 3 + i);
0061       editLineInLineSet(li, p, 3 + i, 0 + i);
0062     }
0063     for (int i = 0; i < 4; ++i)
0064       editLineInLineSet(li, p, 0 + i, 4 + i);
0065   }
0066 }  // namespace
0067
0068 //______________________________________________________________________________
0069 void FWPFCandidateWithHitsProxyBuilder::build(const FWEventItem* iItem,
0070                                               TEveElementList* product,
0071                                               const FWViewContext* vc) {
0072   // init PFCandiate collection
0073   reco::PFCandidateCollection const* candidates = nullptr;
0074   iItem->get(candidates);
0075   if (candidates == nullptr)
0076     return;
0077
0078   Int_t idx = 0;
0079   initPFRecHitsCollections();
0080   for (reco::PFCandidateCollection::const_iterator it = candidates->begin(), itEnd = candidates->end(); it != itEnd;
0081        ++it, ++idx) {
0082     TEveCompound* comp = createCompound();
0083     setupAddElement(comp, product);
0084     // printf("products size %d/%d \n", (int)iItem->size(), product->NumChildren());
0085
0086     const reco::PFCandidate& cand = *it;
0087
0088     // track
0089     {
0090       TEveRecTrack t;
0091       t.fBeta = 1.;
0092       t.fP = TEveVector(cand.px(), cand.py(), cand.pz());
0093       t.fV = TEveVector(cand.vertex().x(), cand.vertex().y(), cand.vertex().z());
0094       t.fSign = cand.charge();
0095       TEveTrack* trk = new TEveTrack(&t, context().getTrackPropagator());
0096       trk->MakeTrack();
0097       fireworks::setTrackTypePF(cand, trk);
0098       setupAddElement(trk, comp);
0099     }
0100     // hits
0101     {
0102       comp->SetMainColor(iItem->defaultDisplayProperties().color());
0103       addHitsForCandidate(cand, comp, vc);
0104     }
0105   }
0106 }
0107
0108 //______________________________________________________________________________
0109 void FWPFCandidateWithHitsProxyBuilder::initPFRecHitsCollections() {
0110   // ref hcal collections
0111   edm::Handle<reco::PFRecHitCollection> handle_hits;
0112
0113   m_collectionHCAL = nullptr;
0114   try {
0115     // edm::InputTag tag("particleFlowRecHitHCAL");
0116     edm::InputTag tag(cname);
0117     item()->getEvent()->getByLabel(tag, handle_hits);
0118     if (handle_hits.isValid()) {
0119       m_collectionHCAL = &*handle_hits;
0120     } else {
0121       fwLog(fwlog::kError) << "FWPFCandidateWithHitsProxyBuilder, item " << item()->name()
0122                            << ": Failed to access collection with name " << cname << "." << std::endl;
0123     }
0124   } catch (...) {
0125     fwLog(fwlog::kError) << "FWPFCandidateWithHitsProxyBuilder, item " << item()->name()
0126                          << ": Failed to access collection with name " << cname << "." << std::endl;
0127   }
0128 }
0129
0130 //______________________________________________________________________________
0131 void FWPFCandidateWithHitsProxyBuilder::viewContextBoxScale(
0132     const float* corners, float scale, bool plotEt, std::vector<float>& scaledCorners, const reco::PFRecHit*) {
0133   static TEveVector vtmp;
0134   vtmp.Set(0.f, 0.f, 0.f);
0135   for (unsigned int i = 0; i < 24; i += 3) {
0136     vtmp[0] += corners[i];
0137     vtmp[1] += corners[i + 1];
0138     vtmp[2] += corners[i + 2];
0139   }
0140   vtmp *= 1.f / 8.f;
0141
0142   if (plotEt) {
0143     scale *= vtmp.Perp() / vtmp.Mag();
0144   }
0145
0146   // Coordinates for a scaled version of the original box
0147   for (unsigned int i = 0; i < 24; i += 3) {
0148     scaledCorners[i] = vtmp[0] + (corners[i] - vtmp[0]) * scale;
0149     scaledCorners[i + 1] = vtmp[1] + (corners[i + 1] - vtmp[1]) * scale;
0150     scaledCorners[i + 2] = vtmp[2] + (corners[i + 2] - vtmp[2]) * scale;
0151   }
0152 }
0153
0154 //______________________________________________________________________________
0155 const reco::PFRecHit* FWPFCandidateWithHitsProxyBuilder::getHitForDetId(unsigned candIdx) {
0156   for (reco::PFRecHitCollection::const_iterator it = m_collectionHCAL->begin(); it != m_collectionHCAL->end(); ++it) {
0157     if (it->detId() == candIdx) {
0158       return &(*it);
0159     }
0160   }
0161   return nullptr;
0162 }
0163
0164 //______________________________________________________________________________
0165 void FWPFCandidateWithHitsProxyBuilder::scaleProduct(TEveElementList* parent,
0166                                                      FWViewType::EType type,
0167                                                      const FWViewContext* vc) {
0168   std::vector<float> scaledCorners(24);
0169
0170   float scale = vc->getEnergyScale()->getScaleFactor3D() / 50;
0171   for (TEveElement::List_i i = parent->BeginChildren(); i != parent->EndChildren(); ++i) {
0172     if ((*i)->NumChildren() > 1) {
0173       TEveElement::List_i xx = (*i)->BeginChildren();
0174       ++xx;
0175       TEveBoxSet* boxset = dynamic_cast<TEveBoxSet*>(*xx);
0176       ++xx;
0177       TEveStraightLineSet* lineset = dynamic_cast<TEveStraightLineSet*>(*xx);
0178       TEveChunkManager::iterator li(lineset->GetLinePlex());
0179       li.next();
0180
0181       TEveChunkManager* plex = boxset->GetPlex();
0182       if (plex->N()) {
0183         for (int atomIdx = 0; atomIdx < plex->Size(); ++atomIdx) {
0184           TEveBoxSet::BFreeBox_t* atom = (TEveBoxSet::BFreeBox_t*)boxset->GetPlex()->Atom(atomIdx);
0185           reco::PFRecHit* hit = (reco::PFRecHit*)boxset->GetUserData(atomIdx);
0186           const float* corners = item()->getGeom()->getCorners(hit->detId());
0187           viewContextBoxScale(corners, hit->energy() * scale, vc->getEnergyScale()->getPlotEt(), scaledCorners, hit);
0188           memcpy(atom->fVertices, &scaledCorners[0], sizeof(atom->fVertices));
0189
0190           editBoxInLineSet(li, &scaledCorners[0]);
0191         }
0192
0193         for (TEveProjectable::ProjList_i p = lineset->BeginProjecteds(); p != lineset->EndProjecteds(); ++p) {
0194           TEveStraightLineSetProjected* projLineSet = (TEveStraightLineSetProjected*)(*p);
0195           projLineSet->UpdateProjection();
0196         }
0197       }
0198     }
0199   }
0200 }
0201 //______________________________________________________________________________
0202 namespace {
0203   TString boxset_tooltip_callback(TEveDigitSet* ds, Int_t idx) {
0204     void* ud = ds->GetUserData(idx);
0205     if (ud) {
0206       reco::PFRecHit* hit = (reco::PFRecHit*)ud;
0207       // printf("idx %d %p hit data %p\n", idx, (void*)hit, ud);
0208       if (hit)
0209         return TString::Format("RecHit %d energy '%f'", idx, hit->energy());
0210       else
0211         return "ERROR";
0212     }
0213     return "NULL";
0214   }
0215 }  // namespace
0216 //______________________________________________________________________________
0217 void FWPFCandidateWithHitsProxyBuilder::addHitsForCandidate(const reco::PFCandidate& cand,
0218                                                             TEveElement* holder,
0219                                                             const FWViewContext* vc) {
0220   reco::PFCandidate::ElementsInBlocks eleInBlocks = cand.elementsInBlocks();
0221
0222   TEveBoxSet* boxset = nullptr;
0223   TEveStraightLineSet* lineset = nullptr;
0224
0225   for (unsigned elIdx = 0; elIdx < eleInBlocks.size(); elIdx++) {
0226     // unsigned ieTrack = 0;
0227     // unsigned ieECAL = 0;
0228     unsigned ieHCAL = 0;
0229
0230     reco::PFBlockRef blockRef = eleInBlocks[elIdx].first;
0231     unsigned indexInBlock = eleInBlocks[elIdx].second;
0232     edm::Ptr<reco::PFBlock> myBlock(blockRef.id(), blockRef.get(), blockRef.key());
0233     /*
0234         if (myBlock->elements()[indexInBlock].type() == 1)
0235         ieTrack = indexInBlock;
0236         if (myBlock->elements()[indexInBlock].type() == 4)
0237         ieECAL = indexInBlock;
0238       */
0239     if (myBlock->elements()[indexInBlock].type() == 5)
0240       ieHCAL = indexInBlock;
0241
0242     std::vector<float> scaledCorners(24);
0243     float scale = vc->getEnergyScale()->getScaleFactor3D() / 50;
0244     if (ieHCAL && m_collectionHCAL) {
0245       reco::PFClusterRef hcalclusterRef = myBlock->elements()[ieHCAL].clusterRef();
0246       edm::Ptr<reco::PFCluster> myCluster(hcalclusterRef.id(), hcalclusterRef.get(), hcalclusterRef.key());
0247       if (myCluster.get()) {
0248         const std::vector<std::pair<DetId, float> >& hitsandfracs = myCluster->hitsAndFractions();
0249
0250         if (!boxset) {
0251           boxset = new TEveBoxSet();
0252           boxset->Reset(TEveBoxSet::kBT_FreeBox, true, hitsandfracs.size());
0253           boxset->SetAntiFlick(true);
0254           boxset->SetAlwaysSecSelect(true);
0255           boxset->SetPickable(true);
0256           boxset->SetTooltipCBFoo(boxset_tooltip_callback);
0257         }
0258
0259         if (!lineset) {
0260           lineset = new TEveStraightLineSet();
0261         }
0262
0263         bool hitsFound = false;
0264         for (int ihandf = 0, lastIdx = (int)(hitsandfracs.size()); ihandf < lastIdx; ihandf++) {
0265           unsigned int hitDetId = hitsandfracs[ihandf].first;
0266           const float* corners = context().getGeom()->getCorners(hitDetId);
0267           const reco::PFRecHit* hit = getHitForDetId(hitDetId);
0268           if (hit) {
0269             viewContextBoxScale(corners, hit->energy() * scale, vc->getEnergyScale()->getPlotEt(), scaledCorners, hit);
0270             boxset->AddBox(&scaledCorners[0]);
0271             // setup last box
0272             boxset->DigitColor(holder->GetMainColor());
0273             boxset->DigitUserData((void*)hit);
0274             addBoxAsLines(lineset, &scaledCorners[0]);
0275             hitsFound = true;
0276           }
0277           /*
0278                // AMT: don't add lines if hit is not found becuse of unconsistency of scaling.
0279                else
0280                {
0281                   addBoxAsLines(lineset, corners);
0282                }
0283                */
0284         }
0285         if (!hitsFound)
0286           fwLog(fwlog::kWarning) << Form(
0287               "Can't find matching hits with for HCAL block %d in %s collection. Number of hits %d.\n",
0288               elIdx,
0289               cname.c_str(),
0290               (int)hitsandfracs.size());
0291
0292       } else {
0293         fwLog(fwlog::kInfo) << "empty cluster \n";
0294       }
0295     }
0296   }  // endloop cand.elementsInBlocks();
0297
0298   if (boxset) {
0299     boxset->RefitPlex();
0300     setupAddElement(boxset, holder);
0301   }
0302
0303   if (lineset) {
0304     setupAddElement(lineset, holder);
0305   }
0306 }
0307
0308 REGISTER_FWPROXYBUILDER(FWPFCandidateWithHitsProxyBuilder,
0309                         reco::PFCandidateCollection,
0310                         "PFCandidatesWithHits",
0311                         FWViewType::kAll3DBits | FWViewType::kAllRPZBits);