File indexing completed on 2024-04-06 12:06:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <memory>
0022
0023
0024
0025 #include <vector>
0026 #include <string>
0027 #include <numeric>
0028
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0031 #include "FWCore/Framework/interface/ConsumesCollector.h"
0032
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/Run.h"
0035 #include "FWCore/Framework/interface/MakerMacros.h"
0036
0037 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0038
0039 #include "FWCore/ServiceRegistry/interface/Service.h"
0040 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0041 #include "TH1F.h"
0042 #include "TH2F.h"
0043
0044 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0045
0046 #include "FWCore/Utilities/interface/InputTag.h"
0047
0048 #include "DataFormats/DetId/interface/DetId.h"
0049 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0050 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0053 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0054
0055 #include "MagneticField/Engine/interface/MagneticField.h"
0056 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0057
0058 #include "FWCore/Framework/interface/ESHandle.h"
0059 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0060 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0061
0062 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0063 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0064
0065
0066
0067
0068
0069 class SeedMultiplicityAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0070 public:
0071 explicit SeedMultiplicityAnalyzer(const edm::ParameterSet&);
0072 ~SeedMultiplicityAnalyzer() override;
0073
0074 class FromTrackRefSeedFilter {
0075 public:
0076 FromTrackRefSeedFilter();
0077 FromTrackRefSeedFilter(edm::ConsumesCollector&& iC, const edm::ParameterSet& iConfig);
0078 const std::string& suffix() const;
0079 void prepareEvent(const edm::Event& iEvent);
0080 bool isSelected(const unsigned int iseed) const;
0081
0082 private:
0083 std::string m_suffix;
0084 bool m_passthrough;
0085 edm::EDGetTokenT<reco::TrackCollection> m_trackcollToken;
0086 edm::EDGetTokenT<reco::TrackRefVector> m_seltrackrefcollToken;
0087 edm::Handle<reco::TrackCollection> m_tracks;
0088 edm::Handle<reco::TrackRefVector> m_seltrackrefs;
0089 };
0090
0091 private:
0092 void analyze(const edm::Event&, const edm::EventSetup&) override;
0093
0094
0095
0096 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> _magFieldToken;
0097 edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> _TTRHBuilderToken;
0098 std::vector<edm::EDGetTokenT<TrajectorySeedCollection>> _seedcollTokens;
0099 std::vector<unsigned int> _seedbins;
0100 std::vector<double> _seedmax;
0101 std::vector<FromTrackRefSeedFilter> _seedfilters;
0102 std::vector<edm::EDGetTokenT<std::map<unsigned int, int>>> _multiplicityMapTokens;
0103 std::vector<std::string> _labels;
0104 std::vector<unsigned int> _selections;
0105 std::vector<unsigned int> _binsmult;
0106 std::vector<unsigned int> _binseta;
0107 std::vector<double> _maxs;
0108 std::vector<TH1F*> _hseedmult;
0109 std::vector<TH1F*> _hseedeta;
0110 std::vector<TH2F*> _hseedphieta;
0111 std::vector<TH1F*> _hpixelrhmult;
0112 std::vector<TH2F*> _hbpixclusleneta;
0113 std::vector<TH2F*> _hfpixclusleneta;
0114 std::vector<TH2F*> _hbpixcluslenangle;
0115 std::vector<TH2F*> _hfpixcluslenangle;
0116 std::vector<std::vector<TH2F*>> _hseedmult2D;
0117 std::vector<std::vector<TH2F*>> _hseedeta2D;
0118 };
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 SeedMultiplicityAnalyzer::SeedMultiplicityAnalyzer(const edm::ParameterSet& iConfig)
0132 : _magFieldToken(esConsumes(edm::ESInputTag{"", iConfig.getParameter<std::string>("TTRHBuilder")})),
0133 _TTRHBuilderToken(esConsumes()),
0134 _seedcollTokens(),
0135 _seedbins(),
0136 _seedmax(),
0137 _seedfilters(),
0138 _multiplicityMapTokens(),
0139 _labels(),
0140 _selections(),
0141 _binsmult(),
0142 _binseta(),
0143 _maxs(),
0144 _hseedmult2D(),
0145 _hseedeta2D() {
0146
0147 usesResource(TFileService::kSharedResource);
0148
0149
0150
0151 std::vector<edm::ParameterSet> seedCollectionConfigs =
0152 iConfig.getParameter<std::vector<edm::ParameterSet>>("seedCollections");
0153
0154 for (std::vector<edm::ParameterSet>::const_iterator scps = seedCollectionConfigs.begin();
0155 scps != seedCollectionConfigs.end();
0156 ++scps) {
0157 _seedcollTokens.push_back(consumes<TrajectorySeedCollection>(scps->getParameter<edm::InputTag>("src")));
0158 _seedbins.push_back(scps->getUntrackedParameter<unsigned int>("nBins", 1000));
0159 _seedmax.push_back(scps->getUntrackedParameter<double>("maxValue", 100000.));
0160
0161 if (scps->exists("trackFilter")) {
0162 _seedfilters.push_back(
0163 FromTrackRefSeedFilter(consumesCollector(), scps->getParameter<edm::ParameterSet>("trackFilter")));
0164 } else {
0165 _seedfilters.push_back(FromTrackRefSeedFilter());
0166 }
0167 }
0168
0169 std::vector<edm::ParameterSet> correlationConfigs =
0170 iConfig.getParameter<std::vector<edm::ParameterSet>>("multiplicityCorrelations");
0171
0172 for (std::vector<edm::ParameterSet>::const_iterator ps = correlationConfigs.begin(); ps != correlationConfigs.end();
0173 ++ps) {
0174 _multiplicityMapTokens.push_back(
0175 consumes<std::map<unsigned int, int>>(ps->getParameter<edm::InputTag>("multiplicityMap")));
0176 _labels.push_back(ps->getParameter<std::string>("detLabel"));
0177 _selections.push_back(ps->getParameter<unsigned int>("detSelection"));
0178 _binsmult.push_back(ps->getParameter<unsigned int>("nBins"));
0179 _binseta.push_back(ps->getParameter<unsigned int>("nBinsEta"));
0180 _maxs.push_back(ps->getParameter<double>("maxValue"));
0181 }
0182
0183 edm::Service<TFileService> tfserv;
0184
0185 std::vector<unsigned int>::const_iterator nseedbins = _seedbins.begin();
0186 std::vector<double>::const_iterator seedmax = _seedmax.begin();
0187 std::vector<FromTrackRefSeedFilter>::const_iterator filter = _seedfilters.begin();
0188
0189 for (std::vector<edm::ParameterSet>::const_iterator scps = seedCollectionConfigs.begin();
0190 scps != seedCollectionConfigs.end();
0191 ++scps, ++nseedbins, ++seedmax, ++filter) {
0192 std::string extendedlabel = std::string(scps->getParameter<edm::InputTag>("src").encode()) + filter->suffix();
0193
0194 std::string hname = extendedlabel + std::string("_mult");
0195 std::string htitle = extendedlabel + std::string(" seed multiplicity");
0196 _hseedmult.push_back(tfserv->make<TH1F>(
0197 hname.c_str(), htitle.c_str(), *nseedbins + 1, 0.5 - *seedmax / (*nseedbins), *seedmax + 0.5));
0198 _hseedmult[_hseedmult.size() - 1]->GetXaxis()->SetTitle("seeds");
0199 _hseedmult[_hseedmult.size() - 1]->GetYaxis()->SetTitle("events");
0200
0201 hname = extendedlabel + std::string("_eta");
0202 htitle = extendedlabel + std::string(" seed pseudorapidity");
0203 _hseedeta.push_back(tfserv->make<TH1F>(hname.c_str(), htitle.c_str(), 80, -4., 4.));
0204 _hseedeta[_hseedeta.size() - 1]->GetXaxis()->SetTitle("#eta");
0205 _hseedeta[_hseedeta.size() - 1]->GetYaxis()->SetTitle("seeds");
0206
0207 hname = extendedlabel + std::string("_phieta");
0208 htitle = extendedlabel + std::string(" seed phi vs pseudorapidity");
0209 _hseedphieta.push_back(tfserv->make<TH2F>(hname.c_str(), htitle.c_str(), 80, -4., 4., 80, -M_PI, M_PI));
0210 _hseedphieta[_hseedphieta.size() - 1]->GetXaxis()->SetTitle("#eta");
0211 _hseedphieta[_hseedphieta.size() - 1]->GetYaxis()->SetTitle("#phi");
0212
0213 _hseedmult2D.push_back(std::vector<TH2F*>());
0214 _hseedeta2D.push_back(std::vector<TH2F*>());
0215
0216 hname = extendedlabel + std::string("_npixelrh");
0217 htitle = extendedlabel + std::string(" seed SiPixelRecHit multiplicity");
0218 _hpixelrhmult.push_back(tfserv->make<TH1F>(hname.c_str(), htitle.c_str(), 5, -.5, 4.5));
0219 _hpixelrhmult[_hpixelrhmult.size() - 1]->GetXaxis()->SetTitle("NRecHits");
0220 _hpixelrhmult[_hpixelrhmult.size() - 1]->GetYaxis()->SetTitle("seeds");
0221
0222 hname = extendedlabel + std::string("_bpixleneta");
0223 htitle = extendedlabel + std::string(" seed BPIX cluster length vs pseudorapidity");
0224 _hbpixclusleneta.push_back(tfserv->make<TH2F>(hname.c_str(), htitle.c_str(), 80, -4., 4., 40, -0.5, 39.5));
0225 _hbpixclusleneta[_hbpixclusleneta.size() - 1]->GetXaxis()->SetTitle("#eta");
0226 _hbpixclusleneta[_hbpixclusleneta.size() - 1]->GetYaxis()->SetTitle("length");
0227
0228 hname = extendedlabel + std::string("_fpixleneta");
0229 htitle = extendedlabel + std::string(" seed FPIX cluster length vs pseudorapidity");
0230 _hfpixclusleneta.push_back(tfserv->make<TH2F>(hname.c_str(), htitle.c_str(), 80, -4., 4., 40, -0.5, 39.5));
0231 _hfpixclusleneta[_hfpixclusleneta.size() - 1]->GetXaxis()->SetTitle("#eta");
0232 _hfpixclusleneta[_hfpixclusleneta.size() - 1]->GetYaxis()->SetTitle("length");
0233
0234 hname = extendedlabel + std::string("_bpixlenangle");
0235 htitle = extendedlabel + std::string(" seed BPIX cluster length vs track projection");
0236 _hbpixcluslenangle.push_back(tfserv->make<TH2F>(hname.c_str(), htitle.c_str(), 200, -1., 1., 40, -0.5, 39.5));
0237 _hbpixcluslenangle[_hbpixcluslenangle.size() - 1]->GetXaxis()->SetTitle("projection");
0238 _hbpixcluslenangle[_hbpixcluslenangle.size() - 1]->GetYaxis()->SetTitle("length");
0239
0240 hname = extendedlabel + std::string("_fpixlenangle");
0241 htitle = extendedlabel + std::string(" seed FPIX cluster length vs track projection");
0242 _hfpixcluslenangle.push_back(tfserv->make<TH2F>(hname.c_str(), htitle.c_str(), 200, -1., 1., 40, -0.5, 39.5));
0243 _hfpixcluslenangle[_hfpixcluslenangle.size() - 1]->GetXaxis()->SetTitle("projection");
0244 _hfpixcluslenangle[_hfpixcluslenangle.size() - 1]->GetYaxis()->SetTitle("length");
0245
0246 for (unsigned int i = 0; i < _multiplicityMapTokens.size(); ++i) {
0247 std::string hname2D = extendedlabel + _labels[i];
0248 hname2D += "_mult";
0249 std::string htitle2D = extendedlabel + " seeds multiplicity";
0250 htitle2D += " vs ";
0251 htitle2D += _labels[i];
0252 htitle2D += " hits";
0253 _hseedmult2D[_hseedmult2D.size() - 1].push_back(tfserv->make<TH2F>(hname2D.c_str(),
0254 htitle2D.c_str(),
0255 _binsmult[i],
0256 0.,
0257 _maxs[i],
0258 *nseedbins + 1,
0259 0.5 - *seedmax / (*nseedbins),
0260 *seedmax + 0.5));
0261 _hseedmult2D[_hseedmult2D.size() - 1][_hseedmult2D[_hseedmult2D.size() - 1].size() - 1]->GetXaxis()->SetTitle(
0262 "hits");
0263 _hseedmult2D[_hseedmult2D.size() - 1][_hseedmult2D[_hseedmult2D.size() - 1].size() - 1]->GetYaxis()->SetTitle(
0264 "seeds");
0265
0266 hname2D = extendedlabel + _labels[i];
0267 hname2D += "_eta";
0268 htitle2D = extendedlabel + " seeds pseudorapidity";
0269 htitle2D += " vs ";
0270 htitle2D += _labels[i];
0271 htitle2D += " hits";
0272 _hseedeta2D[_hseedeta2D.size() - 1].push_back(
0273 tfserv->make<TH2F>(hname2D.c_str(), htitle2D.c_str(), _binseta[i], 0., _maxs[i], 80, -4., 4.));
0274 _hseedeta2D[_hseedeta2D.size() - 1][_hseedeta2D[_hseedeta2D.size() - 1].size() - 1]->GetXaxis()->SetTitle("hits");
0275 _hseedeta2D[_hseedeta2D.size() - 1][_hseedeta2D[_hseedeta2D.size() - 1].size() - 1]->GetYaxis()->SetTitle("#eta");
0276 }
0277 }
0278 }
0279
0280 SeedMultiplicityAnalyzer::~SeedMultiplicityAnalyzer() {
0281
0282
0283 }
0284
0285
0286
0287
0288
0289
0290 void SeedMultiplicityAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0291 using namespace edm;
0292
0293
0294
0295 std::vector<int> tmpmult(_multiplicityMapTokens.size(), -1);
0296 for (unsigned int i = 0; i < _multiplicityMapTokens.size(); ++i) {
0297 Handle<std::map<unsigned int, int>> mults;
0298 iEvent.getByToken(_multiplicityMapTokens[i], mults);
0299
0300
0301
0302 std::map<unsigned int, int>::const_iterator mult = mults->find(_selections[i]);
0303
0304 if (mult != mults->end()) {
0305 tmpmult[i] = mult->second;
0306 } else {
0307 edm::LogWarning("DetSelectionNotFound") << " DetSelection " << _selections[i] << " not found";
0308 }
0309 }
0310
0311
0312
0313
0314 TSCBLBuilderNoMaterial tscblBuilder;
0315
0316 const auto theMF = &iSetup.getData(_magFieldToken);
0317 const auto& theTTRHBuilder = iSetup.getData(_TTRHBuilderToken);
0318
0319
0320
0321
0322 std::vector<TH1F*>::iterator histomult = _hseedmult.begin();
0323 std::vector<std::vector<TH2F*>>::iterator histomult2D = _hseedmult2D.begin();
0324 std::vector<TH1F*>::iterator histoeta = _hseedeta.begin();
0325 std::vector<TH2F*>::iterator histophieta = _hseedphieta.begin();
0326 std::vector<std::vector<TH2F*>>::iterator histoeta2D = _hseedeta2D.begin();
0327 std::vector<TH1F*>::iterator hpixelrhmult = _hpixelrhmult.begin();
0328 std::vector<TH2F*>::iterator histobpixleneta = _hbpixclusleneta.begin();
0329 std::vector<TH2F*>::iterator histofpixleneta = _hfpixclusleneta.begin();
0330 std::vector<TH2F*>::iterator histobpixlenangle = _hbpixcluslenangle.begin();
0331 std::vector<TH2F*>::iterator histofpixlenangle = _hfpixcluslenangle.begin();
0332 std::vector<FromTrackRefSeedFilter>::iterator filter = _seedfilters.begin();
0333
0334
0335
0336 for (std::vector<edm::EDGetTokenT<TrajectorySeedCollection>>::const_iterator coll = _seedcollTokens.begin();
0337 coll != _seedcollTokens.end() && histomult != _hseedmult.end() && histomult2D != _hseedmult2D.end() &&
0338 histoeta != _hseedeta.end() && histoeta2D != _hseedeta2D.end() && histophieta != _hseedphieta.end() &&
0339 hpixelrhmult != _hpixelrhmult.end() && histobpixleneta != _hbpixclusleneta.end() &&
0340 histofpixleneta != _hfpixclusleneta.end() && histobpixlenangle != _hbpixcluslenangle.end() &&
0341 histofpixlenangle != _hfpixcluslenangle.end();
0342 ++coll,
0343 ++histomult,
0344 ++histomult2D,
0345 ++histoeta,
0346 ++histophieta,
0347 ++histoeta2D,
0348 ++hpixelrhmult,
0349 ++histobpixleneta,
0350 ++histofpixleneta,
0351 ++histobpixlenangle,
0352 ++histofpixlenangle,
0353 ++filter) {
0354 filter->prepareEvent(iEvent);
0355
0356 Handle<TrajectorySeedCollection> seeds;
0357 iEvent.getByToken(*coll, seeds);
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369 unsigned int nseeds = 0;
0370 unsigned int iseed = 0;
0371 for (TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed, ++iseed) {
0372 if (filter->isSelected(iseed)) {
0373 ++nseeds;
0374
0375 TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder.build(&*(seed->recHits().end() - 1));
0376 TrajectoryStateOnSurface state =
0377 trajectoryStateTransform::transientState(seed->startingState(), recHit->surface(), theMF);
0378
0379
0380 double eta = state.globalMomentum().eta();
0381 double phi = state.globalMomentum().phi();
0382
0383 (*histoeta)->Fill(eta);
0384 (*histophieta)->Fill(eta, phi);
0385
0386 for (unsigned int i = 0; i < _multiplicityMapTokens.size(); ++i) {
0387 if (tmpmult[i] >= 0)
0388 (*histoeta2D)[i]->Fill(tmpmult[i], eta);
0389 }
0390
0391 int npixelrh = 0;
0392 for (auto const& hit : seed->recHits()) {
0393 const SiPixelRecHit* sphit = dynamic_cast<const SiPixelRecHit*>(&hit);
0394 if (sphit) {
0395 ++npixelrh;
0396
0397 TransientTrackingRecHit::RecHitPointer ttrhit = theTTRHBuilder.build(&hit);
0398 TrajectoryStateOnSurface tsos =
0399 trajectoryStateTransform::transientState(seed->startingState(), ttrhit->surface(), theMF);
0400
0401 if (sphit->geographicalId().det() == DetId::Tracker &&
0402 sphit->geographicalId().subdetId() == PixelSubdetector::PixelBarrel) {
0403 (*histobpixleneta)->Fill(eta, sphit->cluster()->sizeY());
0404 if (tsos.isValid()) {
0405
0406 double normdx = tsos.localMomentum().x() / sqrt(tsos.localMomentum().x() * tsos.localMomentum().x() +
0407 tsos.localMomentum().z() * tsos.localMomentum().z());
0408 (*histobpixlenangle)->Fill(normdx, sphit->cluster()->sizeY());
0409 }
0410 } else if (sphit->geographicalId().det() == DetId::Tracker &&
0411 sphit->geographicalId().subdetId() == PixelSubdetector::PixelEndcap) {
0412 (*histofpixleneta)->Fill(eta, sphit->cluster()->sizeX());
0413 if (tsos.isValid()) {
0414
0415 double normdy = tsos.localMomentum().y() / sqrt(tsos.localMomentum().y() * tsos.localMomentum().y() +
0416 tsos.localMomentum().z() * tsos.localMomentum().z());
0417 (*histofpixlenangle)->Fill(normdy, sphit->cluster()->sizeX());
0418 }
0419 } else {
0420 edm::LogError("InconsistentSiPixelRecHit")
0421 << "SiPixelRecHit with a non-pixel DetId " << sphit->geographicalId().rawId();
0422 }
0423 }
0424 }
0425 (*hpixelrhmult)->Fill(npixelrh);
0426 }
0427 }
0428 (*histomult)->Fill(nseeds);
0429
0430 for (unsigned int i = 0; i < _multiplicityMapTokens.size(); ++i) {
0431 if (tmpmult[i] >= 0)
0432 (*histomult2D)[i]->Fill(tmpmult[i], nseeds);
0433 }
0434 }
0435 }
0436
0437 SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::FromTrackRefSeedFilter()
0438 : m_suffix(""), m_passthrough(true), m_trackcollToken(), m_seltrackrefcollToken(), m_tracks(), m_seltrackrefs() {}
0439
0440 SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::FromTrackRefSeedFilter(edm::ConsumesCollector&& iC,
0441 const edm::ParameterSet& iConfig)
0442 : m_suffix(iConfig.getParameter<std::string>("suffix")),
0443 m_passthrough(false),
0444 m_trackcollToken(iC.consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trkCollection"))),
0445 m_seltrackrefcollToken(
0446 iC.consumes<reco::TrackRefVector>(iConfig.getParameter<edm::InputTag>("selRefTrkCollection"))),
0447 m_tracks(),
0448 m_seltrackrefs() {}
0449
0450 const std::string& SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::suffix() const { return m_suffix; }
0451
0452 void SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::prepareEvent(const edm::Event& iEvent) {
0453 if (!m_passthrough) {
0454 iEvent.getByToken(m_trackcollToken, m_tracks);
0455 iEvent.getByToken(m_seltrackrefcollToken, m_seltrackrefs);
0456 }
0457
0458 return;
0459 }
0460
0461 bool SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::isSelected(const unsigned int iseed) const {
0462 if (m_passthrough) {
0463 return true;
0464 } else {
0465
0466 const reco::TrackRef trkref(m_tracks, iseed);
0467
0468
0469 for (reco::TrackRefVector::const_iterator seltrkref = m_seltrackrefs->begin(); seltrkref != m_seltrackrefs->end();
0470 ++seltrkref) {
0471 if (trkref == *seltrkref)
0472 return true;
0473 }
0474 }
0475 return false;
0476 }
0477
0478
0479 DEFINE_FWK_MODULE(SeedMultiplicityAnalyzer);