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   //loop on tracks
0125   for (std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk;
0126        ++ittrk) {
0127     //loop on tracking rechits
0128     LogDebug("AlignmentPrescaler") << "Loop on hits of track #" << (ittrk - Tracks->begin()) << std::endl;
0129     int ntakenhits = 0;
0130     bool firstTakenHit = false;
0131
0132     for (auto const& hit : ittrk->recHits()) {
0133       if (!hit->isValid()) {
0134         continue;
0135       }
0136       uint32_t tmpdetid = hit->geographicalId().rawId();
0137       tpresc_->GetEntryWithIndex(tmpdetid);
0138
0139       //-------------
0140       //decide whether to take this hit or not
0141       bool takeit = false;
0142       int subdetId = hit->geographicalId().subdetId();
0143
0144       //check first if the cluster is also in the overlap asso map
0145       bool isOverlapHit = false;
0146       //  bool first=true;
0147       //ugly...
0148       const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0149       const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0150       const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0151
0152       AlignmentClusterFlag tmpflag(hit->geographicalId());
0153       int stripType = 0;
0154       if (subdetId > 2) {  // SST case
0155         const std::type_info& type = typeid(*hit);
0156         if (type == typeid(SiStripRecHit1D))
0157           stripType = 1;
0158         else if (type == typeid(SiStripRecHit2D))
0159           stripType = 2;
0160         else
0161           stripType = 3;
0162
0163         if (stripType == 1) {
0164           //      const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0165
0166           if (stripHit1D != nullptr) {
0167             SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
0168             tmpflag = InValMap[stripclust];
0169             tmpflag.SetDetId(hit->geographicalId());
0170             if (tmpflag.isOverlap())
0171               isOverlapHit = true;
0172             LogDebug("AlignmentPrescaler")
0173                 << "~*~*~* Prescale (1D) for module " << tmpflag.detId().rawId() << "("
0174                 << InValMap[stripclust].detId().rawId() << ") is " << hitPrescFactor_ << std::flush;
0175             if (tmpflag.isOverlap())
0176               LogDebug("AlignmentPrescaler") << " (it is Overlap)";
0177           }  //end if striphit1D!=0
0178         } else if (stripType == 2) {
0179           //const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0180           if (stripHit2D != nullptr) {
0181             SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
0182             tmpflag = InValMap[stripclust];
0183             tmpflag.SetDetId(hit->geographicalId());
0184             if (tmpflag.isOverlap())
0185               isOverlapHit = true;
0186             LogDebug("AlignmentPrescaler")
0187                 << "~*~*~* Prescale (2D) for module " << tmpflag.detId().rawId() << "("
0188                 << InValMap[stripclust].detId().rawId() << ") is " << hitPrescFactor_ << std::flush;
0189             if (tmpflag.isOverlap())
0190               LogDebug("AlignmentPrescaler") << " (it is Overlap)" << endl;
0191           }  //end if striphit2D!=0
0192         }
0193       }  //end if is a strip hit
0194       else {
0195         //  const SiPixelRecHit*   pixelhit= dynamic_cast<const SiPixelRecHit*>(hit);
0196         if (pixelhit != nullptr) {
0197           SiPixelClusterRefNew pixclust(pixelhit->cluster());
0198           tmpflag = InValMap[pixclust];
0199           tmpflag.SetDetId(hit->geographicalId());
0200           if (tmpflag.isOverlap())
0201             isOverlapHit = true;
0202         }
0203       }  //end else is a pixel hit
0204       //      tmpflag.SetDetId(hit->geographicalId());
0205
0206       if (isOverlapHit) {
0207         LogDebug("AlignmentPrescaler") << "  DetId=" << tmpdetid << " is Overlap! " << flush;
0208         takeit = (float(myrand_->Rndm()) <= overlapPrescFactor_);
0209       }
0210       if (!takeit) {
0211         float rr = float(myrand_->Rndm());
0212         takeit = (rr <= hitPrescFactor_);
0213       }
0214       if (takeit) {  //HIT TAKEN !
0215         LogDebug("AlignmentPrescaler") << "  DetId=" << tmpdetid << " taken!" << flush;
0216         tmpflag.SetTakenFlag();
0217
0218         if (subdetId > 2) {
0219           if (stripType == 1) {
0220             SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
0221             InValMap[stripclust] = tmpflag;  //.SetTakenFlag();
0222           } else if (stripType == 2) {
0223             SiStripRecHit1D::ClusterRef stripclust(stripHit2D->cluster());
0224             InValMap[stripclust] = tmpflag;  //.SetTakenFlag();
0225           } else
0226             std::cout << "Unknown type of strip hit" << std::endl;
0227         } else {
0228           SiPixelClusterRefNew pixclust(pixelhit->cluster());
0229           InValMap[pixclust] = tmpflag;  //.SetTakenFlag();
0230         }
0231
0232         if (!firstTakenHit) {
0233           firstTakenHit = true;
0234           LogDebug("AlignmentPrescaler") << "Index of the track iterator is " << ittrk - Tracks->begin();
0235         }
0236         ntakenhits++;
0237       }  //end if take this hit
0238     }  //end loop on RecHits
0239     trackflags[ittrk - Tracks->begin()] = ntakenhits;
0240   }  //end loop on tracks
0241
0242   //save the asso map, tracks...
0243   // prepare output
0244   auto OutVM = std::make_unique<AliClusterValueMap>();
0245   *OutVM = InValMap;
0246
0247   iEvent.put(std::move(OutVM));
0248
0249   auto trkVM = std::make_unique<AliTrackTakenClusterValueMap>();
0250   AliTrackTakenClusterValueMap::Filler trkmapfiller(*trkVM);
0251   trkmapfiller.insert(Tracks, trackflags.begin(), trackflags.end());
0252   trkmapfiller.fill();
0253   iEvent.put(std::move(trkVM));
0254
0255 }  //end produce
0256
0257 int AlignmentPrescaler::layerFromId(const DetId& id, const TrackerTopology* tTopo) const {
0258   if (uint32_t(id.subdetId()) == PixelSubdetector::PixelBarrel) {
0259     return tTopo->pxbLayer(id);
0260   } else if (uint32_t(id.subdetId()) == PixelSubdetector::PixelEndcap) {
0261     return tTopo->pxfDisk(id) + (3 * (tTopo->pxfSide(id) - 1));
0262   } else if (id.subdetId() == StripSubdetector::TIB) {
0263     return tTopo->tibLayer(id);
0264   } else if (id.subdetId() == StripSubdetector::TOB) {
0265     return tTopo->tobLayer(id);
0266   } else if (id.subdetId() == StripSubdetector::TEC) {
0267     return tTopo->tecWheel(id) + (9 * (tTopo->pxfSide(id) - 1));
0268   } else if (id.subdetId() == StripSubdetector::TID) {
0269     return tTopo->tidWheel(id) + (3 * (tTopo->tidSide(id) - 1));
0270   }
0271   return -1;
0272 }  //end layerfromId
0273
0274 void AlignmentPrescaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0275   edm::ParameterSetDescription desc;
0276   desc.setComment("Prescale Tracker Alignment hits for alignment procedures");
0277   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0278   desc.add<edm::InputTag>("assomap", edm::InputTag("OverlapAssoMap"));
0279   desc.add<std::string>("PrescFileName", "PrescaleFactors.root");
0280   desc.add<std::string>("PrescTreeName", "AlignmentHitMap");
0281   descriptions.addWithDefaultLabel(desc);
0282 }
0283
0284 // ========= MODULE DEF ==============
0285 #include "FWCore/PluginManager/interface/ModuleDef.h"
0286 #include "FWCore/Framework/interface/MakerMacros.h"
0287 DEFINE_FWK_MODULE(AlignmentPrescaler);