TrackerAlignment/plugins/TkAlCaOverlapTagger.cc

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/EventPrincipal.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011
0012 #include "DataFormats/Common/interface/View.h"
0013 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0014 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0015 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0016
0017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0020 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0021
0022 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0023 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0024 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0025 #include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h"
0026 #include "Utilities/General/interface/ClassName.h"
0027
0028 #include "DataFormats/Alignment/interface/AlignmentClusterFlag.h"
0029 #include "DataFormats/Alignment/interface/AliClusterValueMap.h"
0030
0031 class TkAlCaOverlapTagger : public edm::stream::EDProducer<> {
0032 public:
0033   TkAlCaOverlapTagger(const edm::ParameterSet& iConfig);
0034   ~TkAlCaOverlapTagger() override;
0035   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0036   static void fillDescriptions(edm::ConfigurationDescriptions&);
0037
0038 private:
0039   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0040   edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackToken_;
0041   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> siPixelClustersToken_;
0042   edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> siStripClustersToken_;
0043   edm::InputTag src_;
0044   edm::InputTag srcClust_;
0045   bool rejectBadMods_;
0046   std::vector<unsigned int> BadModsList_;
0047
0048   int layerFromId(const DetId& id, const TrackerTopology* tTopo) const;
0049 };
0050
0051 TkAlCaOverlapTagger::TkAlCaOverlapTagger(const edm::ParameterSet& iConfig)
0052     : topoToken_(esConsumes()),
0053       trajTrackToken_(consumes((iConfig.getParameter<edm::InputTag>("src")))),
0054       siPixelClustersToken_(consumes((iConfig.getParameter<edm::InputTag>("Clustersrc")))),
0055       siStripClustersToken_(consumes((iConfig.getParameter<edm::InputTag>("Clustersrc")))),
0056       rejectBadMods_(iConfig.getParameter<bool>("rejectBadMods")),
0057       BadModsList_(iConfig.getParameter<std::vector<uint32_t>>("BadMods")) {
0058   produces<AliClusterValueMap>();  //produces the ValueMap (VM) to be stored in the Event at the end
0059 }
0060
0061 TkAlCaOverlapTagger::~TkAlCaOverlapTagger() = default;
0062
0063 void TkAlCaOverlapTagger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0064   //Retrieve tracker topology from geometry
0065   const TrackerTopology* tTopo = &iSetup.getData(topoToken_);
0066
0067   edm::Handle<TrajTrackAssociationCollection> assoMap;
0068   iEvent.getByToken(trajTrackToken_, assoMap);
0069   LogDebug("TkAlCaOverlapTagger") << "\n\n############################\n###  Starting a new TkAlCaOverlapTagger - Ev "
0070                                   << iEvent.id().run() << ", " << iEvent.id().event();
0071
0072   AlignmentClusterFlag iniflag;
0073   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> pixelclusters;
0074   iEvent.getByToken(siPixelClustersToken_, pixelclusters);  //same label as tracks
0075   std::vector<AlignmentClusterFlag> pixelvalues(pixelclusters->dataSize(),
0076                                                 iniflag);  //vector where to store value to be fileld in the VM
0077
0078   edm::Handle<edmNew::DetSetVector<SiStripCluster>> stripclusters;
0079   iEvent.getByToken(siStripClustersToken_, stripclusters);  //same label as tracks
0080   std::vector<AlignmentClusterFlag> stripvalues(stripclusters->dataSize(),
0081                                                 iniflag);  //vector where to store value to be fileld in the VM
0082
0083   //start doing the thing!
0084
0085   //loop over trajectories
0086   for (TrajTrackAssociationCollection::const_iterator itass = assoMap->begin(); itass != assoMap->end(); ++itass) {
0087     int nOverlaps = 0;
0088     const edm::Ref<std::vector<Trajectory>> traj = itass->key;  //trajectory in the collection
0089     const Trajectory* myTrajectory = &(*traj);
0090     std::vector<TrajectoryMeasurement> tmColl = myTrajectory->measurements();
0091
0092     const reco::TrackRef tkref = itass->val;  //associated track track in the collection
0093     // const Track * trk = &(*tkref);
0094     int hitcnt = -1;
0095
0096     //loop over traj meas
0097     const TrajectoryMeasurement* previousTM(nullptr);
0098     DetId previousId(0);
0099
0100     for (std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas != tmColl.end();
0101          ++itTrajMeas) {
0102       hitcnt++;
0103
0104       if (previousTM != nullptr) {
0105         LogDebug("TkAlCaOverlapTagger") << "Checking TrajMeas (" << hitcnt + 1 << "):";
0106         if (!previousTM->recHit()->isValid()) {
0107           LogDebug("TkAlCaOverlapTagger") << "Previous RecHit invalid !";
0108           continue;
0109         } else {
0110           LogDebug("TkAlCaOverlapTagger") << "\nDetId: " << std::flush << previousTM->recHit()->geographicalId().rawId()
0111                                           << "\t Local x of hit: " << previousTM->recHit()->localPosition().x();
0112         }
0113       } else {
0114         LogDebug("TkAlCaOverlapTagger") << "This is the first Traj Meas of the Trajectory! The Trajectory contains "
0115                                         << tmColl.size() << " TrajMeas";
0116       }
0117
0118       TrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit();
0119       const TrackingRecHit* hit = &(*hitpointer);
0120       if (!hit->isValid())
0121         continue;
0122
0123       DetId detid = hit->geographicalId();
0124       int layer(layerFromId(detid, tTopo));  //layer 1-4=TIB, layer 5-10=TOB
0125       int subDet = detid.subdetId();
0126
0127       if ((previousTM != nullptr) && (layer != -1)) {
0128         for (std::vector<TrajectoryMeasurement>::const_iterator itmCompare = itTrajMeas - 1;
0129              itmCompare >= tmColl.begin() && itmCompare > itTrajMeas - 4;
0130              --itmCompare) {
0131           DetId compareId = itmCompare->recHit()->geographicalId();
0132           if (subDet != compareId.subdetId() || layer != layerFromId(compareId, tTopo))
0133             break;
0134           if (!itmCompare->recHit()->isValid())
0135             continue;
0136           if ((subDet <= 2) ||
0137               (subDet > 2 && SiStripDetId(detid).stereo() ==
0138                                  SiStripDetId(compareId).stereo())) {  //if either pixel or strip stereo module
0139
0140             //
0141             //WOW, we have an overlap!!!!!!
0142             //
0143             AlignmentClusterFlag hitflag(hit->geographicalId());
0144             hitflag.SetOverlapFlag();
0145             LogDebug("TkAlCaOverlapTagger") << "Overlap found in SubDet " << subDet << "!!!" << std::flush;
0146
0147             bool hitInStrip = (subDet == SiStripDetId::TIB) || (subDet == SiStripDetId::TID) ||
0148                               (subDet == SiStripDetId::TOB) || (subDet == SiStripDetId::TEC);
0149             if (hitInStrip) {
0150               LogDebug("TkAlCaOverlapTagger") << "  TypeId of the RecHit: " << className(*hit);
0151               // const std::type_info &type = typeid(*hit);
0152               const SiStripRecHit2D* transstriphit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0153               const SiStripRecHit1D* transstriphit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0154
0155               //   if (type == typeid(SiStripRecHit1D)) {
0156               if (transstriphit1D != nullptr) {
0157                 //  const SiStripRecHit1D* striphit=dynamic_cast<const  SiStripRecHit1D*>(hit);
0158                 const SiStripRecHit1D* striphit = transstriphit1D;
0159                 if (striphit != nullptr) {
0160                   SiStripRecHit1D::ClusterRef stripclust(striphit->cluster());
0161
0162                   if (stripclust.id() ==
0163                       stripclusters
0164                           .id()) {  //ensure that the stripclust is really present in the original cluster collection!!!
0165                     stripvalues[stripclust.key()] = hitflag;
0166                   } else {
0167                     edm::LogError("TkAlCaOverlapTagger")
0168                         << "ERROR in <TkAlCaOverlapTagger::produce>: ProdId of Strip clusters mismatched: "
0169                         << stripclust.id() << " vs " << stripclusters.id();
0170                   }
0171                 } else {
0172                   edm::LogError("TkAlCaOverlapTagger") << "ERROR in <TkAlCaOverlapTagger::produce>: Dynamic cast of "
0173                                                           "Strip RecHit failed!   TypeId of the RecHit: "
0174                                                        << className(*hit);
0175                 }
0176               }  //end if sistriprechit1D
0177               else if (transstriphit2D != nullptr) {
0178                 //else if (type == typeid(SiStripRecHit2D)) {
0179                 //      const SiStripRecHit2D* striphit=dynamic_cast<const  SiStripRecHit2D*>(hit);
0180                 const SiStripRecHit2D* striphit = transstriphit2D;
0181                 if (striphit != nullptr) {
0182                   SiStripRecHit2D::ClusterRef stripclust(striphit->cluster());
0183
0184                   if (stripclust.id() ==
0185                       stripclusters
0186                           .id()) {  //ensure that the stripclust is really present in the original cluster collection!!!
0187                     stripvalues[stripclust.key()] = hitflag;
0188
0189                     LogDebug("TkAlCaOverlapTagger")
0190                         << ">>> Storing in the ValueMap a StripClusterRef with Cluster.Key: " << stripclust.key()
0191                         << " (" << striphit->cluster().key() << "), Cluster.Id: " << stripclust.id() << "  (DetId is "
0192                         << hit->geographicalId().rawId() << ")";
0193                   } else {
0194                     edm::LogError("TkAlCaOverlapTagger")
0195                         << "ERROR in <TkAlCaOverlapTagger::produce>: ProdId of Strip clusters mismatched: "
0196                         << stripclust.id() << " vs " << stripclusters.id();
0197                   }
0198
0199                   LogDebug("TkAlCaOverlapTagger") << "Cluster baricentre: " << stripclust->barycenter();
0200                 } else {
0201                   edm::LogError("TkAlCaOverlapTagger") << "ERROR in <TkAlCaOverlapTagger::produce>: Dynamic cast of "
0202                                                           "Strip RecHit failed!   TypeId of the RecHit: "
0203                                                        << className(*hit);
0204                 }
0205               }  //end if Sistriprechit2D
0206               else {
0207                 edm::LogError("TkAlCaOverlapTagger") << "ERROR in <TkAlCaOverlapTagger::produce>: Impossible to "
0208                                                         "determine the type of SiStripRecHit.  TypeId of the RecHit: "
0209                                                      << className(*hit);
0210               }
0211
0212             }       //end if hit in Strips
0213             else {  //pixel hit
0214               const SiPixelRecHit* transpixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0215               if (transpixelhit != nullptr) {
0216                 const SiPixelRecHit* pixelhit = transpixelhit;
0217                 SiPixelClusterRefNew pixclust(pixelhit->cluster());
0218
0219                 if (pixclust.id() == pixelclusters.id()) {
0220                   pixelvalues[pixclust.key()] = hitflag;
0221                   LogDebug("TkAlCaOverlapTagger")
0222                       << ">>> Storing in the ValueMap a PixelClusterRef with ProdID: " << pixclust.id()
0223                       << "  (DetId is " << hit->geographicalId().rawId() << ")";  //"  and  a Val with ID: "<<flag.id();
0224                 } else {
0225                   edm::LogError("TkAlCaOverlapTagger")
0226                       << "ERROR in <TkAlCaOverlapTagger::produce>: ProdId of Pixel clusters mismatched: "
0227                       << pixclust.id() << " vs " << pixelclusters.id();
0228                 }
0229               } else {
0230                 edm::LogError("TkAlCaOverlapTagger") << "ERROR in <TkAlCaOverlapTagger::produce>: Dynamic cast of "
0231                                                         "Pixel RecHit failed!   TypeId of the RecHit: "
0232                                                      << className(*hit);
0233               }
0234             }  //end 'else' it is a pixel hit
0235
0236             nOverlaps++;
0237             break;
0238           }
0239         }  //end second loop on TM
0240       }    //end if a previous TM exists
0241
0242       previousTM = &(*itTrajMeas);
0243       previousId = detid;
0244     }  //end loop over traj meas
0245     LogDebug("TkAlCaOverlapTagger") << "Found " << nOverlaps << " overlaps in this trajectory";
0246   }  //end loop over trajectories
0247
0248   // prepare output
0249   auto hitvalmap = std::make_unique<AliClusterValueMap>();
0250   AliClusterValueMap::Filler mapfiller(*hitvalmap);
0251
0252   edm::TestHandle<std::vector<AlignmentClusterFlag>> fakePixelHandle(&pixelvalues, pixelclusters.id());
0253   mapfiller.insert(fakePixelHandle, pixelvalues.begin(), pixelvalues.end());
0254
0255   edm::TestHandle<std::vector<AlignmentClusterFlag>> fakeStripHandle(&stripvalues, stripclusters.id());
0256   mapfiller.insert(fakeStripHandle, stripvalues.begin(), stripvalues.end());
0257   mapfiller.fill();
0258
0259   // iEvent.put(std::move(stripmap));
0260   iEvent.put(std::move(hitvalmap));
0261 }  //end  TkAlCaOverlapTagger::produce
0262 int TkAlCaOverlapTagger::layerFromId(const DetId& id, const TrackerTopology* tTopo) const {
0263   if (uint32_t(id.subdetId()) == PixelSubdetector::PixelBarrel) {
0264     return tTopo->pxbLayer(id);
0265   } else if (uint32_t(id.subdetId()) == PixelSubdetector::PixelEndcap) {
0266     return tTopo->pxfDisk(id) + (3 * (tTopo->pxfSide(id) - 1));
0267   } else if (id.subdetId() == StripSubdetector::TIB) {
0268     return tTopo->tibLayer(id);
0269   } else if (id.subdetId() == StripSubdetector::TOB) {
0270     return tTopo->tobLayer(id);
0271   } else if (id.subdetId() == StripSubdetector::TEC) {
0272     return tTopo->tecWheel(id) + (9 * (tTopo->tecSide(id) - 1));
0273   } else if (id.subdetId() == StripSubdetector::TID) {
0274     return tTopo->tidWheel(id) + (3 * (tTopo->tidSide(id) - 1));
0275   }
0276   return -1;
0277
0278 }  //end layerfromId
0279
0280 void TkAlCaOverlapTagger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0281   edm::ParameterSetDescription desc;
0282   desc.setComment("Tagger of overlaps for tracker alignment");
0283   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0284   desc.add<edm::InputTag>("Clustersrc", edm::InputTag("ALCARECOTkAlCosmicsCTF0T"));
0285   desc.add<bool>("rejectBadMods", false);
0286   desc.add<std::vector<uint32_t>>("BadMods", {});
0287   descriptions.addWithDefaultLabel(desc);
0288 }
0289
0290 // ========= MODULE DEF ==============
0291 #include "FWCore/PluginManager/interface/ModuleDef.h"
0292 #include "FWCore/Framework/interface/MakerMacros.h"
0293
0294 DEFINE_FWK_MODULE(TkAlCaOverlapTagger);