TrackerAlignment/plugins/AlignmentPrescaler.cc

0001 #include "Alignment/CommonAlignment/interface/Alignable.h"
0002 #include "Alignment/CommonAlignment/interface/Utilities.h"
0003 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0004 #include "AnalysisDataFormats/TrackInfo/interface/TrackInfo.h"
0005 #include "DataFormats/Alignment/interface/AliClusterValueMap.h"
0006 #include "DataFormats/Alignment/interface/AlignmentClusterFlag.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0012 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventPrincipal.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/one/EDProducer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026 #include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h"
0027 #include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h"
0028 #include "Utilities/General/interface/ClassName.h"
0029
0030 #include "TFile.h"
0031 #include "TTree.h"
0032 #include "TRandom3.h"
0033 #include "TH1F.h"
0034 #include <string>
0035
0036 class TrackerTopology;
0037
0038 class AlignmentPrescaler : public edm::one::EDProducer<> {
0039 public:
0040   AlignmentPrescaler(const edm::ParameterSet& iConfig);
0041   ~AlignmentPrescaler() override;
0042   void beginJob() override;
0043   void endJob() override;
0044   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0045   static void fillDescriptions(edm::ConfigurationDescriptions&);
0046
0047 private:
0048   // tokens
0049   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;   //tracks in input
0050   edm::EDGetTokenT<AliClusterValueMap> aliClusterToken_;  //Hit-quality association map
0051
0052   std::string prescfilename_;  //name of the file containing the TTree with the prescaling factors
0053   std::string presctreename_;  //name of the  TTree with the prescaling factors
0054
0055   TFile* fpresc_;
0056   TTree* tpresc_;
0057   TRandom3* myrand_;
0058
0059   int layerFromId(const DetId& id, const TrackerTopology* tTopo) const;
0060
0061   unsigned int detid_;
0062   float hitPrescFactor_, overlapPrescFactor_;
0063   int totnhitspxl_;
0064 };
0065
0066 using namespace std;
0067
0068 AlignmentPrescaler::AlignmentPrescaler(const edm::ParameterSet& iConfig)
0069     : tracksToken_(consumes(iConfig.getParameter<edm::InputTag>("src"))),
0070       aliClusterToken_(consumes(iConfig.getParameter<edm::InputTag>("assomap"))),
0071       prescfilename_(iConfig.getParameter<std::string>("PrescFileName")),
0072       presctreename_(iConfig.getParameter<std::string>("PrescTreeName")) {
0073   // issue the produce<>
0074   produces<AliClusterValueMap>();
0075   produces<AliTrackTakenClusterValueMap>();
0076 }
0077
0078 AlignmentPrescaler::~AlignmentPrescaler() = default;
0079
0080 void AlignmentPrescaler::beginJob() {
0081   //
0082   edm::LogPrint("AlignmentPrescaler") << "in AlignmentPrescaler::beginJob" << std::flush;
0083   fpresc_ = new TFile(prescfilename_.c_str(), "READ");
0084   tpresc_ = (TTree*)fpresc_->Get(presctreename_.c_str());
0085   tpresc_->BuildIndex("DetId");
0086   tpresc_->SetBranchStatus("*", false);
0087   tpresc_->SetBranchStatus("DetId", true);
0088   tpresc_->SetBranchStatus("PrescaleFactor", true);
0089   tpresc_->SetBranchStatus("PrescaleFactorOverlap", true);
0090   edm::LogPrint("AlignmentPrescaler") << " Branches activated " << std::flush;
0091   detid_ = 0;
0092   hitPrescFactor_ = 99.0;
0093   overlapPrescFactor_ = 88.0;
0094
0095   tpresc_->SetBranchAddress("DetId", &detid_);
0096   tpresc_->SetBranchAddress("PrescaleFactor", &hitPrescFactor_);
0097   tpresc_->SetBranchAddress("PrescaleFactorOverlap", &overlapPrescFactor_);
0098   edm::LogPrint("AlignmentPrescaler") << " addressed " << std::flush;
0099   myrand_ = new TRandom3();
0100   //   myrand_->SetSeed();
0101   edm::LogPrint("AlignmentPrescaler") << " ok ";
0102 }
0103
0104 void AlignmentPrescaler::endJob() {
0105   delete tpresc_;
0106   fpresc_->Close();
0107   delete fpresc_;
0108   delete myrand_;
0109 }
0110
0111 void AlignmentPrescaler::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0112   LogDebug("AlignmentPrescaler") << "\n\n#################\n### Starting the AlignmentPrescaler::produce ; Event: "
0113                                  << iEvent.id().run() << ", " << iEvent.id().event() << std::endl;
0114
0115   edm::Handle<reco::TrackCollection> Tracks;
0116   iEvent.getByToken(tracksToken_, Tracks);
0117
0118   //take  HitAssomap
0119   AliClusterValueMap InValMap = iEvent.get(aliClusterToken_);
0120
0121   //prepare the output of the ValueMap flagging tracks
0122   std::vector<int> trackflags(Tracks->size(), 0);
0123
0124   int npxlhits = 0;
0125
0126   //loop on tracks
0127   for (std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk;
0128        ++ittrk) {
0129     //loop on tracking rechits
0130     LogDebug("AlignmentPrescaler") << "Loop on hits of track #" << (ittrk - Tracks->begin()) << std::endl;
0131     int nhit = 0;
0132     int ntakenhits = 0;
0133     bool firstTakenHit = false;
0134
0135     for (auto const& hit : ittrk->recHits()) {
0136       if (!hit->isValid()) {
0137         nhit++;
0138         continue;
0139       }
0140       uint32_t tmpdetid = hit->geographicalId().rawId();
0141       tpresc_->GetEntryWithIndex(tmpdetid);
0142
0143       //-------------
0144       //decide whether to take this hit or not
0145       bool takeit = false;
0146       int subdetId = hit->geographicalId().subdetId();
0147
0148       //check first if the cluster is also in the overlap asso map
0149       bool isOverlapHit = false;
0150       //  bool first=true;
0151       //ugly...
0152       const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0153       const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0154       const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0155
0156       AlignmentClusterFlag tmpflag(hit->geographicalId());
0157       int stripType = 0;
0158       if (subdetId > 2) {  // SST case
0159         const std::type_info& type = typeid(*hit);
0160         if (type == typeid(SiStripRecHit1D))
0161           stripType = 1;
0162         else if (type == typeid(SiStripRecHit2D))
0163           stripType = 2;
0164         else
0165           stripType = 3;
0166
0167         if (stripType == 1) {
0168           //      const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0169
0170           if (stripHit1D != nullptr) {
0171             SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
0172             tmpflag = InValMap[stripclust];
0173             tmpflag.SetDetId(hit->geographicalId());
0174             if (tmpflag.isOverlap())
0175               isOverlapHit = true;
0176             LogDebug("AlignmentPrescaler")
0177                 << "~*~*~* Prescale (1D) for module " << tmpflag.detId().rawId() << "("
0178                 << InValMap[stripclust].detId().rawId() << ") is " << hitPrescFactor_ << std::flush;
0179             if (tmpflag.isOverlap())
0180               LogDebug("AlignmentPrescaler") << " (it is Overlap)";
0181           }  //end if striphit1D!=0
0182         } else if (stripType == 2) {
0183           //const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0184           if (stripHit2D != nullptr) {
0185             SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
0186             tmpflag = InValMap[stripclust];
0187             tmpflag.SetDetId(hit->geographicalId());
0188             if (tmpflag.isOverlap())
0189               isOverlapHit = true;
0190             LogDebug("AlignmentPrescaler")
0191                 << "~*~*~* Prescale (2D) for module " << tmpflag.detId().rawId() << "("
0192                 << InValMap[stripclust].detId().rawId() << ") is " << hitPrescFactor_ << std::flush;
0193             if (tmpflag.isOverlap())
0194               LogDebug("AlignmentPrescaler") << " (it is Overlap)" << endl;
0195           }  //end if striphit2D!=0
0196         }
0197       }  //end if is a strip hit
0198       else {
0199         //  const SiPixelRecHit*   pixelhit= dynamic_cast<const SiPixelRecHit*>(hit);
0200         if (pixelhit != nullptr) {
0201           npxlhits++;
0202           SiPixelClusterRefNew pixclust(pixelhit->cluster());
0203           tmpflag = InValMap[pixclust];
0204           tmpflag.SetDetId(hit->geographicalId());
0205           if (tmpflag.isOverlap())
0206             isOverlapHit = true;
0207         }
0208       }  //end else is a pixel hit
0209       //      tmpflag.SetDetId(hit->geographicalId());
0210
0211       if (isOverlapHit) {
0212         LogDebug("AlignmentPrescaler") << "  DetId=" << tmpdetid << " is Overlap! " << flush;
0213         takeit = (float(myrand_->Rndm()) <= overlapPrescFactor_);
0214       }
0215       if (!takeit) {
0216         float rr = float(myrand_->Rndm());
0217         takeit = (rr <= hitPrescFactor_);
0218       }
0219       if (takeit) {  //HIT TAKEN !
0220         LogDebug("AlignmentPrescaler") << "  DetId=" << tmpdetid << " taken!" << flush;
0221         tmpflag.SetTakenFlag();
0222
0223         if (subdetId > 2) {
0224           if (stripType == 1) {
0225             SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
0226             InValMap[stripclust] = tmpflag;  //.SetTakenFlag();
0227           } else if (stripType == 2) {
0228             SiStripRecHit1D::ClusterRef stripclust(stripHit2D->cluster());
0229             InValMap[stripclust] = tmpflag;  //.SetTakenFlag();
0230           } else
0231             std::cout << "Unknown type of strip hit" << std::endl;
0232         } else {
0233           SiPixelClusterRefNew pixclust(pixelhit->cluster());
0234           InValMap[pixclust] = tmpflag;  //.SetTakenFlag();
0235         }
0236
0237         if (!firstTakenHit) {
0238           firstTakenHit = true;
0239           LogDebug("AlignmentPrescaler") << "Index of the track iterator is " << ittrk - Tracks->begin();
0240         }
0241         ntakenhits++;
0242       }  //end if take this hit
0243       nhit++;
0244     }  //end loop on RecHits
0245     trackflags[ittrk - Tracks->begin()] = ntakenhits;
0246   }  //end loop on tracks
0247
0248   //save the asso map, tracks...
0249   // prepare output
0250   auto OutVM = std::make_unique<AliClusterValueMap>();
0251   *OutVM = InValMap;
0252
0253   iEvent.put(std::move(OutVM));
0254
0255   auto trkVM = std::make_unique<AliTrackTakenClusterValueMap>();
0256   AliTrackTakenClusterValueMap::Filler trkmapfiller(*trkVM);
0257   trkmapfiller.insert(Tracks, trackflags.begin(), trackflags.end());
0258   trkmapfiller.fill();
0259   iEvent.put(std::move(trkVM));
0260
0261 }  //end produce
0262
0263 int AlignmentPrescaler::layerFromId(const DetId& id, const TrackerTopology* tTopo) const {
0264   if (uint32_t(id.subdetId()) == PixelSubdetector::PixelBarrel) {
0265     return tTopo->pxbLayer(id);
0266   } else if (uint32_t(id.subdetId()) == PixelSubdetector::PixelEndcap) {
0267     return tTopo->pxfDisk(id) + (3 * (tTopo->pxfSide(id) - 1));
0268   } else if (id.subdetId() == StripSubdetector::TIB) {
0269     return tTopo->tibLayer(id);
0270   } else if (id.subdetId() == StripSubdetector::TOB) {
0271     return tTopo->tobLayer(id);
0272   } else if (id.subdetId() == StripSubdetector::TEC) {
0273     return tTopo->tecWheel(id) + (9 * (tTopo->pxfSide(id) - 1));
0274   } else if (id.subdetId() == StripSubdetector::TID) {
0275     return tTopo->tidWheel(id) + (3 * (tTopo->tidSide(id) - 1));
0276   }
0277   return -1;
0278 }  //end layerfromId
0279
0280 void AlignmentPrescaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0281   edm::ParameterSetDescription desc;
0282   desc.setComment("Prescale Tracker Alignment hits for alignment procedures");
0283   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0284   desc.add<edm::InputTag>("assomap", edm::InputTag("OverlapAssoMap"));
0285   desc.add<std::string>("PrescFileName", "PrescaleFactors.root");
0286   desc.add<std::string>("PrescTreeName", "AlignmentHitMap");
0287   descriptions.addWithDefaultLabel(desc);
0288 }
0289
0290 // ========= MODULE DEF ==============
0291 #include "FWCore/PluginManager/interface/ModuleDef.h"
0292 #include "FWCore/Framework/interface/MakerMacros.h"
0293 DEFINE_FWK_MODULE(AlignmentPrescaler);