File indexing completed on 2021-09-16 03:23:16
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
0049 edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0050 edm::EDGetTokenT<AliClusterValueMap> aliClusterToken_;
0051
0052 std::string prescfilename_;
0053 std::string presctreename_;
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
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
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
0119 AliClusterValueMap InValMap = iEvent.get(aliClusterToken_);
0120
0121
0122 std::vector<int> trackflags(Tracks->size(), 0);
0123
0124 int npxlhits = 0;
0125
0126
0127 for (std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk;
0128 ++ittrk) {
0129
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
0145 bool takeit = false;
0146 int subdetId = hit->geographicalId().subdetId();
0147
0148
0149 bool isOverlapHit = false;
0150
0151
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) {
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
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 }
0182 } else if (stripType == 2) {
0183
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 }
0196 }
0197 }
0198 else {
0199
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 }
0209
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) {
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;
0227 } else if (stripType == 2) {
0228 SiStripRecHit1D::ClusterRef stripclust(stripHit2D->cluster());
0229 InValMap[stripclust] = tmpflag;
0230 } else
0231 std::cout << "Unknown type of strip hit" << std::endl;
0232 } else {
0233 SiPixelClusterRefNew pixclust(pixelhit->cluster());
0234 InValMap[pixclust] = tmpflag;
0235 }
0236
0237 if (!firstTakenHit) {
0238 firstTakenHit = true;
0239 LogDebug("AlignmentPrescaler") << "Index of the track iterator is " << ittrk - Tracks->begin();
0240 }
0241 ntakenhits++;
0242 }
0243 nhit++;
0244 }
0245 trackflags[ittrk - Tracks->begin()] = ntakenhits;
0246 }
0247
0248
0249
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 }
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 }
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
0291 #include "FWCore/PluginManager/interface/ModuleDef.h"
0292 #include "FWCore/Framework/interface/MakerMacros.h"
0293 DEFINE_FWK_MODULE(AlignmentPrescaler);