File indexing completed on 2021-08-23 03:25:28
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0010 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0012 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0013 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
0014 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0015 #include "DataFormats/GeometrySurface/interface/TkRotation.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0021 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0022 #include "DataFormats/ParticleFlowReco/interface/PFSimParticle.h"
0023 #include "DataFormats/ParticleFlowReco/interface/PFSimParticleFwd.h"
0024 #include "DataFormats/TrackReco/interface/Track.h"
0025 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0026 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/stream/EDProducer.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0034 #include "FWCore/Utilities/interface/EDGetToken.h"
0035 #include "FWCore/Utilities/interface/Exception.h"
0036 #include "FastSimulation/Event/interface/FSimEvent.h"
0037 #include "FastSimulation/Event/interface/FSimTrack.h"
0038 #include "FastSimulation/Event/interface/FSimVertex.h"
0039 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0045 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
0046 #include "RecoParticleFlow/PFProducer/interface/PFBlockAlgo.h"
0047 #include "RecoParticleFlow/PFTracking/interface/PFGeometry.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0049 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0050 #include "SimDataFormats/Track/interface/SimTrack.h"
0051 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0052 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0053 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0054 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0055 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0056 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0057 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0058 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0059
0060 #include <memory>
0061 #include <set>
0062 #include <sstream>
0063 #include <string>
0064
0065 class PFSimParticleProducer : public edm::stream::EDProducer<> {
0066 public:
0067 explicit PFSimParticleProducer(const edm::ParameterSet&);
0068
0069 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0070
0071 void produce(edm::Event&, const edm::EventSetup&) override;
0072
0073 typedef edm::Handle<reco::PFRecTrackCollection> TrackHandle;
0074 void getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID);
0075
0076 private:
0077
0078 edm::InputTag inputTagSim_;
0079 edm::EDGetTokenT<std::vector<SimTrack> > tokenSim_;
0080 edm::EDGetTokenT<std::vector<SimVertex> > tokenSimVertices_;
0081
0082
0083
0084 bool mctruthMatchingInfo_;
0085 edm::InputTag inputTagFastSimProducer_;
0086 edm::EDGetTokenT<edm::PCaloHitContainer> tokenFastSimProducer_;
0087
0088
0089 edm::InputTag inputTagRecTracks_;
0090 edm::EDGetTokenT<reco::PFRecTrackCollection> tokenRecTracks_;
0091 edm::InputTag inputTagEcalRecHitsEB_;
0092 edm::EDGetTokenT<EcalRecHitCollection> tokenEcalRecHitsEB_;
0093 edm::InputTag inputTagEcalRecHitsEE_;
0094 edm::EDGetTokenT<EcalRecHitCollection> tokenEcalRecHitsEE_;
0095
0096 edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;
0097
0098
0099
0100 edm::ParameterSet particleFilter_;
0101 std::unique_ptr<FSimEvent> mySimEvent;
0102
0103
0104
0105
0106 bool processParticles_;
0107
0108
0109 bool verbose_;
0110 };
0111
0112 #include "FWCore/Framework/interface/MakerMacros.h"
0113 DEFINE_FWK_MODULE(PFSimParticleProducer);
0114
0115 void PFSimParticleProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0116
0117 edm::ParameterSetDescription desc;
0118 desc.addUntracked<edm::InputTag>("fastSimProducer", edm::InputTag("fastSimProducer", "EcalHitsEB"));
0119 desc.addUntracked<bool>("MCTruthMatchingInfo", false);
0120 desc.add<edm::InputTag>("RecTracks", edm::InputTag("trackerDrivenElectronSeeds"));
0121 desc.add<std::string>("Fitter", "KFFittingSmoother");
0122 desc.add<edm::InputTag>("ecalRecHitsEE", edm::InputTag("caloRecHits", "EcalRecHitsEE"));
0123 desc.add<edm::InputTag>("ecalRecHitsEB", edm::InputTag("caloRecHits", "EcalRecHitsEB"));
0124 desc.addUntracked<bool>("process_RecTracks", false);
0125 {
0126 edm::ParameterSetDescription psd0;
0127 psd0.setUnknown();
0128 desc.add<edm::ParameterSetDescription>("ParticleFilter", psd0);
0129 }
0130 desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0131 desc.addUntracked<bool>("process_Particles", true);
0132 desc.add<std::string>("Propagator", "PropagatorWithMaterial");
0133 desc.add<edm::InputTag>("sim", edm::InputTag("g4SimHits"));
0134 desc.addUntracked<bool>("verbose", false);
0135 descriptions.add("particleFlowSimParticle", desc);
0136 }
0137
0138 using namespace std;
0139 using namespace edm;
0140
0141 PFSimParticleProducer::PFSimParticleProducer(const edm::ParameterSet& iConfig) {
0142 processParticles_ = iConfig.getUntrackedParameter<bool>("process_Particles", true);
0143
0144 inputTagSim_ = iConfig.getParameter<InputTag>("sim");
0145 tokenSim_ = consumes<std::vector<SimTrack> >(inputTagSim_);
0146 tokenSimVertices_ = consumes<std::vector<SimVertex> >(inputTagSim_);
0147
0148
0149
0150
0151 inputTagFastSimProducer_ = iConfig.getUntrackedParameter<InputTag>("fastSimProducer");
0152 tokenFastSimProducer_ = consumes<edm::PCaloHitContainer>(inputTagFastSimProducer_);
0153 mctruthMatchingInfo_ = iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo", false);
0154
0155
0156 inputTagRecTracks_ = iConfig.getParameter<InputTag>("RecTracks");
0157 tokenRecTracks_ = consumes<reco::PFRecTrackCollection>(inputTagRecTracks_);
0158
0159 inputTagEcalRecHitsEB_ = iConfig.getParameter<InputTag>("ecalRecHitsEB");
0160 tokenEcalRecHitsEB_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEB_);
0161 inputTagEcalRecHitsEE_ = iConfig.getParameter<InputTag>("ecalRecHitsEE");
0162 tokenEcalRecHitsEE_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEE_);
0163
0164 pdtToken_ = esConsumes();
0165
0166 verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0167
0168
0169 produces<reco::PFSimParticleCollection>();
0170
0171 particleFilter_ = iConfig.getParameter<ParameterSet>("ParticleFilter");
0172
0173 mySimEvent = std::make_unique<FSimEvent>(particleFilter_);
0174 }
0175
0176 void PFSimParticleProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0177
0178 mySimEvent->initializePdt(&iSetup.getData(pdtToken_));
0179
0180 LogDebug("PFSimParticleProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0181 << endl;
0182
0183
0184
0185
0186
0187 std::vector<unsigned> recTrackSimID;
0188
0189
0190
0191
0192
0193 typedef std::pair<double, unsigned> hitSimID;
0194 std::vector<std::list<hitSimID> > caloHitsEBID(62000);
0195 std::vector<double> caloHitsEBTotE(62000, 0.0);
0196
0197 if (mctruthMatchingInfo_) {
0198
0199 auto pcalohits = iEvent.getHandle(tokenFastSimProducer_);
0200
0201 if (!pcalohits) {
0202 ostringstream err;
0203 err << "could not find pcaloHit "
0204 << "fastSimProducer:EcalHitsEB";
0205 LogError("PFSimParticleProducer") << err.str() << endl;
0206
0207 throw cms::Exception("MissingProduct", err.str());
0208 } else {
0209 assert(pcalohits.isValid());
0210
0211 edm::PCaloHitContainer::const_iterator it = pcalohits.product()->begin();
0212 edm::PCaloHitContainer::const_iterator itend = pcalohits.product()->end();
0213
0214
0215 for (; it != itend; ++it) {
0216 EBDetId detid(it->id());
0217
0218 if (it->energy() > 0.0) {
0219 std::pair<double, unsigned> phitsimid = make_pair(it->energy(), it->geantTrackId());
0220 caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
0221 caloHitsEBTotE[detid.hashedIndex()] += it->energy();
0222 }
0223
0224 }
0225 }
0226
0227
0228
0229 Handle<reco::PFRecTrackCollection> recTracks;
0230 try {
0231 LogDebug("PFSimParticleProducer") << "getting PFRecTracks" << endl;
0232 iEvent.getByToken(tokenRecTracks_, recTracks);
0233
0234 } catch (cms::Exception& err) {
0235 LogError("PFSimParticleProducer") << err << " cannot get collection "
0236 << "particleFlowBlock"
0237 << ":"
0238 << "" << endl;
0239 }
0240
0241
0242
0243 getSimIDs(recTracks, recTrackSimID);
0244
0245 }
0246
0247
0248 if (processParticles_) {
0249 auto pOutputPFSimParticleCollection = std::make_unique<reco::PFSimParticleCollection>();
0250
0251 Handle<vector<SimTrack> > simTracks;
0252 bool found = iEvent.getByToken(tokenSim_, simTracks);
0253 if (!found) {
0254 ostringstream err;
0255 err << "cannot find sim tracks: " << inputTagSim_;
0256 LogError("PFSimParticleProducer") << err.str() << endl;
0257
0258 throw cms::Exception("MissingProduct", err.str());
0259 }
0260
0261 Handle<vector<SimVertex> > simVertices;
0262 found = iEvent.getByToken(tokenSimVertices_, simVertices);
0263 if (!found) {
0264 LogError("PFSimParticleProducer") << "cannot find sim vertices: " << inputTagSim_ << endl;
0265 return;
0266 }
0267
0268 mySimEvent->fill(*simTracks, *simVertices);
0269
0270 if (verbose_)
0271 mySimEvent->print();
0272
0273 for (unsigned i = 0; i < mySimEvent->nTracks(); i++) {
0274 const FSimTrack& fst = mySimEvent->track(i);
0275
0276 int motherId = -1;
0277 if (!fst.noMother())
0278 motherId = fst.mother().id();
0279
0280
0281
0282
0283
0284 unsigned recTrackID = 99999;
0285 vector<unsigned> recHitContrib;
0286 vector<double> recHitContribFrac;
0287
0288 if (mctruthMatchingInfo_) {
0289
0290 for (unsigned lo = 0; lo < recTrackSimID.size(); lo++) {
0291 if (i == recTrackSimID[lo]) {
0292 recTrackID = lo;
0293 }
0294 }
0295
0296
0297 edm::Handle<EcalRecHitCollection> rhcHandle;
0298 bool found = iEvent.getByToken(tokenEcalRecHitsEB_, rhcHandle);
0299 if (!found) {
0300 ostringstream err;
0301 err << "could not find rechits " << inputTagEcalRecHitsEB_;
0302 LogError("PFSimParticleProducer") << err.str() << endl;
0303
0304 throw cms::Exception("MissingProduct", err.str());
0305 } else {
0306 assert(rhcHandle.isValid());
0307
0308 EBRecHitCollection::const_iterator it_rh = rhcHandle.product()->begin();
0309 EBRecHitCollection::const_iterator itend_rh = rhcHandle.product()->end();
0310
0311 for (; it_rh != itend_rh; ++it_rh) {
0312 unsigned rhit_hi = EBDetId(it_rh->id()).hashedIndex();
0313 EBDetId detid(it_rh->id());
0314
0315 auto it_phit = caloHitsEBID[rhit_hi].begin();
0316 auto itend_phit = caloHitsEBID[rhit_hi].end();
0317 for (; it_phit != itend_phit; ++it_phit) {
0318 if (i == it_phit->second) {
0319
0320
0321 bool alreadyin = false;
0322 for (unsigned ihit = 0; ihit < recHitContrib.size(); ++ihit)
0323 if (detid.rawId() == recHitContrib[ihit])
0324 alreadyin = true;
0325
0326 if (!alreadyin) {
0327 double pcalofraction = 0.0;
0328 if (caloHitsEBTotE[rhit_hi] != 0.0)
0329 pcalofraction = (it_phit->first / caloHitsEBTotE[rhit_hi]) * 100.0;
0330
0331
0332 recHitContrib.push_back(it_rh->id());
0333 recHitContribFrac.push_back(pcalofraction);
0334 }
0335 }
0336 }
0337
0338 }
0339
0340 }
0341
0342 }
0343
0344 reco::PFSimParticle particle(
0345 fst.charge(), fst.type(), fst.id(), motherId, fst.daughters(), recTrackID, recHitContrib, recHitContribFrac);
0346
0347 const FSimVertex& originVtx = fst.vertex();
0348
0349 math::XYZPoint posOrig(originVtx.position().x(), originVtx.position().y(), originVtx.position().z());
0350
0351 math::XYZTLorentzVector momOrig(
0352 fst.momentum().px(), fst.momentum().py(), fst.momentum().pz(), fst.momentum().e());
0353 reco::PFTrajectoryPoint pointOrig(-1, reco::PFTrajectoryPoint::ClosestApproach, posOrig, momOrig);
0354
0355
0356 particle.addPoint(pointOrig);
0357
0358 if (!fst.noEndVertex()) {
0359 const FSimVertex& endVtx = fst.endVertex();
0360
0361 math::XYZPoint posEnd(endVtx.position().x(), endVtx.position().y(), endVtx.position().z());
0362
0363 math::XYZTLorentzVector momEnd;
0364
0365 reco::PFTrajectoryPoint pointEnd(-1, reco::PFTrajectoryPoint::BeamPipeOrEndVertex, posEnd, momEnd);
0366
0367 particle.addPoint(pointEnd);
0368 } else {
0369 reco::PFTrajectoryPoint dummy;
0370 particle.addPoint(dummy);
0371 }
0372
0373 if (fst.onLayer1()) {
0374 const RawParticle& rp = fst.layer1Entrance();
0375
0376 math::XYZPoint posLayer1(rp.x(), rp.y(), rp.z());
0377 math::XYZTLorentzVector momLayer1(rp.px(), rp.py(), rp.pz(), rp.e());
0378 reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1, posLayer1, momLayer1);
0379
0380 particle.addPoint(layer1Pt);
0381
0382
0383 } else {
0384 reco::PFTrajectoryPoint dummy;
0385 particle.addPoint(dummy);
0386 }
0387
0388 if (fst.onLayer2()) {
0389 const RawParticle& rp = fst.layer2Entrance();
0390
0391 math::XYZPoint posLayer2(rp.x(), rp.y(), rp.z());
0392 math::XYZTLorentzVector momLayer2(rp.px(), rp.py(), rp.pz(), rp.e());
0393 reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2, posLayer2, momLayer2);
0394
0395 particle.addPoint(layer2Pt);
0396
0397
0398 } else {
0399 reco::PFTrajectoryPoint dummy;
0400 particle.addPoint(dummy);
0401 }
0402
0403 if (fst.onEcal()) {
0404 const RawParticle& rp = fst.ecalEntrance();
0405
0406 math::XYZPoint posECAL(rp.x(), rp.y(), rp.z());
0407 math::XYZTLorentzVector momECAL(rp.px(), rp.py(), rp.pz(), rp.e());
0408 reco::PFTrajectoryPoint ecalPt(-1, reco::PFTrajectoryPoint::ECALEntrance, posECAL, momECAL);
0409
0410 particle.addPoint(ecalPt);
0411
0412
0413 } else {
0414 reco::PFTrajectoryPoint dummy;
0415 particle.addPoint(dummy);
0416 }
0417
0418
0419 reco::PFTrajectoryPoint dummy;
0420 particle.addPoint(dummy);
0421
0422 if (fst.onHcal()) {
0423 const RawParticle& rpin = fst.hcalEntrance();
0424
0425 math::XYZPoint posHCALin(rpin.x(), rpin.y(), rpin.z());
0426 math::XYZTLorentzVector momHCALin(rpin.px(), rpin.py(), rpin.pz(), rpin.e());
0427 reco::PFTrajectoryPoint hcalPtin(-1, reco::PFTrajectoryPoint::HCALEntrance, posHCALin, momHCALin);
0428
0429 particle.addPoint(hcalPtin);
0430
0431 const RawParticle& rpout = fst.hcalExit();
0432
0433 math::XYZPoint posHCALout(rpout.x(), rpout.y(), rpout.z());
0434 math::XYZTLorentzVector momHCALout(rpout.px(), rpout.py(), rpout.pz(), rpout.e());
0435 reco::PFTrajectoryPoint hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit, posHCALout, momHCALout);
0436
0437 particle.addPoint(hcalPtout);
0438
0439 const RawParticle& rpho = fst.hoEntrance();
0440
0441 math::XYZPoint posHOEntrance(rpho.x(), rpho.y(), rpho.z());
0442 math::XYZTLorentzVector momHOEntrance(rpho.px(), rpho.py(), rpho.pz(), rpho.e());
0443 reco::PFTrajectoryPoint hoPtin(-1, reco::PFTrajectoryPoint::HOLayer, posHOEntrance, momHOEntrance);
0444
0445 particle.addPoint(hoPtin);
0446
0447 } else {
0448 reco::PFTrajectoryPoint dummy;
0449 particle.addPoint(dummy);
0450 }
0451
0452 pOutputPFSimParticleCollection->push_back(particle);
0453 }
0454
0455 iEvent.put(std::move(pOutputPFSimParticleCollection));
0456 }
0457
0458 LogDebug("PFSimParticleProducer") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
0459 }
0460
0461 void PFSimParticleProducer::getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID) {
0462 if (trackh.isValid()) {
0463 for (unsigned i = 0; i < trackh->size(); i++) {
0464 const reco::PFRecTrackRef ref(trackh, i);
0465
0466 for (auto const& hit : ref->trackRef()->recHits()) {
0467 if (hit->isValid()) {
0468 auto rechit = dynamic_cast<const FastTrackerRecHit*>(hit);
0469
0470 for (unsigned int st_index = 0; st_index < rechit->nSimTrackIds(); ++st_index) {
0471 recTrackSimID.push_back(rechit->simTrackId(st_index));
0472 }
0473 break;
0474 }
0475 }
0476 }
0477 }
0478 }