File indexing completed on 2024-09-07 04:37:52
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 desc.add<edm::ParameterSetDescription>("ParticleFilter", {});
0126 desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0127 desc.addUntracked<bool>("process_Particles", true);
0128 desc.add<std::string>("Propagator", "PropagatorWithMaterial");
0129 desc.add<edm::InputTag>("sim", edm::InputTag("g4SimHits"));
0130 desc.addUntracked<bool>("verbose", false);
0131 descriptions.add("particleFlowSimParticle", desc);
0132 }
0133
0134 using namespace std;
0135 using namespace edm;
0136
0137 PFSimParticleProducer::PFSimParticleProducer(const edm::ParameterSet& iConfig) {
0138 processParticles_ = iConfig.getUntrackedParameter<bool>("process_Particles", true);
0139
0140 inputTagSim_ = iConfig.getParameter<InputTag>("sim");
0141 tokenSim_ = consumes<std::vector<SimTrack> >(inputTagSim_);
0142 tokenSimVertices_ = consumes<std::vector<SimVertex> >(inputTagSim_);
0143
0144
0145
0146
0147 inputTagFastSimProducer_ = iConfig.getUntrackedParameter<InputTag>("fastSimProducer");
0148 tokenFastSimProducer_ = consumes<edm::PCaloHitContainer>(inputTagFastSimProducer_);
0149 mctruthMatchingInfo_ = iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo", false);
0150
0151
0152 inputTagRecTracks_ = iConfig.getParameter<InputTag>("RecTracks");
0153 tokenRecTracks_ = consumes<reco::PFRecTrackCollection>(inputTagRecTracks_);
0154
0155 inputTagEcalRecHitsEB_ = iConfig.getParameter<InputTag>("ecalRecHitsEB");
0156 tokenEcalRecHitsEB_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEB_);
0157 inputTagEcalRecHitsEE_ = iConfig.getParameter<InputTag>("ecalRecHitsEE");
0158 tokenEcalRecHitsEE_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEE_);
0159
0160 pdtToken_ = esConsumes();
0161
0162 verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0163
0164
0165 produces<reco::PFSimParticleCollection>();
0166
0167 particleFilter_ = iConfig.getParameter<ParameterSet>("ParticleFilter");
0168
0169 mySimEvent = std::make_unique<FSimEvent>(particleFilter_);
0170 }
0171
0172 void PFSimParticleProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0173
0174 mySimEvent->initializePdt(&iSetup.getData(pdtToken_));
0175
0176 LogDebug("PFSimParticleProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0177 << endl;
0178
0179
0180
0181
0182
0183 std::vector<unsigned> recTrackSimID;
0184
0185
0186
0187
0188
0189 typedef std::pair<double, unsigned> hitSimID;
0190 std::vector<std::list<hitSimID> > caloHitsEBID(62000);
0191 std::vector<double> caloHitsEBTotE(62000, 0.0);
0192
0193 if (mctruthMatchingInfo_) {
0194
0195 auto pcalohits = iEvent.getHandle(tokenFastSimProducer_);
0196
0197 if (!pcalohits) {
0198 ostringstream err;
0199 err << "could not find pcaloHit "
0200 << "fastSimProducer:EcalHitsEB";
0201 LogError("PFSimParticleProducer") << err.str() << endl;
0202
0203 throw cms::Exception("MissingProduct", err.str());
0204 } else {
0205 assert(pcalohits.isValid());
0206
0207 edm::PCaloHitContainer::const_iterator it = pcalohits.product()->begin();
0208 edm::PCaloHitContainer::const_iterator itend = pcalohits.product()->end();
0209
0210
0211 for (; it != itend; ++it) {
0212 EBDetId detid(it->id());
0213
0214 if (it->energy() > 0.0) {
0215 std::pair<double, unsigned> phitsimid = make_pair(it->energy(), it->geantTrackId());
0216 caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
0217 caloHitsEBTotE[detid.hashedIndex()] += it->energy();
0218 }
0219
0220 }
0221 }
0222
0223
0224
0225 Handle<reco::PFRecTrackCollection> recTracks;
0226 try {
0227 LogDebug("PFSimParticleProducer") << "getting PFRecTracks" << endl;
0228 iEvent.getByToken(tokenRecTracks_, recTracks);
0229
0230 } catch (cms::Exception& err) {
0231 LogError("PFSimParticleProducer") << err << " cannot get collection "
0232 << "particleFlowBlock"
0233 << ":"
0234 << "" << endl;
0235 }
0236
0237
0238
0239 getSimIDs(recTracks, recTrackSimID);
0240
0241 }
0242
0243
0244 if (processParticles_) {
0245 auto pOutputPFSimParticleCollection = std::make_unique<reco::PFSimParticleCollection>();
0246
0247 Handle<vector<SimTrack> > simTracks;
0248 bool found = iEvent.getByToken(tokenSim_, simTracks);
0249 if (!found) {
0250 ostringstream err;
0251 err << "cannot find sim tracks: " << inputTagSim_;
0252 LogError("PFSimParticleProducer") << err.str() << endl;
0253
0254 throw cms::Exception("MissingProduct", err.str());
0255 }
0256
0257 Handle<vector<SimVertex> > simVertices;
0258 found = iEvent.getByToken(tokenSimVertices_, simVertices);
0259 if (!found) {
0260 LogError("PFSimParticleProducer") << "cannot find sim vertices: " << inputTagSim_ << endl;
0261 return;
0262 }
0263
0264 mySimEvent->fill(*simTracks, *simVertices);
0265
0266 if (verbose_)
0267 mySimEvent->print();
0268
0269 for (unsigned i = 0; i < mySimEvent->nTracks(); i++) {
0270 const FSimTrack& fst = mySimEvent->track(i);
0271
0272 int motherId = -1;
0273 if (!fst.noMother())
0274 motherId = fst.mother().id();
0275
0276
0277
0278
0279
0280 unsigned recTrackID = 99999;
0281 vector<unsigned> recHitContrib;
0282 vector<double> recHitContribFrac;
0283
0284 if (mctruthMatchingInfo_) {
0285
0286 for (unsigned lo = 0; lo < recTrackSimID.size(); lo++) {
0287 if (i == recTrackSimID[lo]) {
0288 recTrackID = lo;
0289 }
0290 }
0291
0292
0293 edm::Handle<EcalRecHitCollection> rhcHandle;
0294 bool found = iEvent.getByToken(tokenEcalRecHitsEB_, rhcHandle);
0295 if (!found) {
0296 ostringstream err;
0297 err << "could not find rechits " << inputTagEcalRecHitsEB_;
0298 LogError("PFSimParticleProducer") << err.str() << endl;
0299
0300 throw cms::Exception("MissingProduct", err.str());
0301 } else {
0302 assert(rhcHandle.isValid());
0303
0304 EBRecHitCollection::const_iterator it_rh = rhcHandle.product()->begin();
0305 EBRecHitCollection::const_iterator itend_rh = rhcHandle.product()->end();
0306
0307 for (; it_rh != itend_rh; ++it_rh) {
0308 unsigned rhit_hi = EBDetId(it_rh->id()).hashedIndex();
0309 EBDetId detid(it_rh->id());
0310
0311 auto it_phit = caloHitsEBID[rhit_hi].begin();
0312 auto itend_phit = caloHitsEBID[rhit_hi].end();
0313 for (; it_phit != itend_phit; ++it_phit) {
0314 if (i == it_phit->second) {
0315
0316
0317 bool alreadyin = false;
0318 for (unsigned ihit = 0; ihit < recHitContrib.size(); ++ihit)
0319 if (detid.rawId() == recHitContrib[ihit])
0320 alreadyin = true;
0321
0322 if (!alreadyin) {
0323 double pcalofraction = 0.0;
0324 if (caloHitsEBTotE[rhit_hi] != 0.0)
0325 pcalofraction = (it_phit->first / caloHitsEBTotE[rhit_hi]) * 100.0;
0326
0327
0328 recHitContrib.push_back(it_rh->id());
0329 recHitContribFrac.push_back(pcalofraction);
0330 }
0331 }
0332 }
0333
0334 }
0335
0336 }
0337
0338 }
0339
0340 reco::PFSimParticle particle(
0341 fst.charge(), fst.type(), fst.id(), motherId, fst.daughters(), recTrackID, recHitContrib, recHitContribFrac);
0342
0343 const FSimVertex& originVtx = fst.vertex();
0344
0345 math::XYZPoint posOrig(originVtx.position().x(), originVtx.position().y(), originVtx.position().z());
0346
0347 math::XYZTLorentzVector momOrig(
0348 fst.momentum().px(), fst.momentum().py(), fst.momentum().pz(), fst.momentum().e());
0349 reco::PFTrajectoryPoint pointOrig(-1, reco::PFTrajectoryPoint::ClosestApproach, posOrig, momOrig);
0350
0351
0352 particle.addPoint(pointOrig);
0353
0354 if (!fst.noEndVertex()) {
0355 const FSimVertex& endVtx = fst.endVertex();
0356
0357 math::XYZPoint posEnd(endVtx.position().x(), endVtx.position().y(), endVtx.position().z());
0358
0359 math::XYZTLorentzVector momEnd;
0360
0361 reco::PFTrajectoryPoint pointEnd(-1, reco::PFTrajectoryPoint::BeamPipeOrEndVertex, posEnd, momEnd);
0362
0363 particle.addPoint(pointEnd);
0364 } else {
0365 reco::PFTrajectoryPoint dummy;
0366 particle.addPoint(dummy);
0367 }
0368
0369 if (fst.onLayer1()) {
0370 const RawParticle& rp = fst.layer1Entrance();
0371
0372 math::XYZPoint posLayer1(rp.x(), rp.y(), rp.z());
0373 math::XYZTLorentzVector momLayer1(rp.px(), rp.py(), rp.pz(), rp.e());
0374 reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1, posLayer1, momLayer1);
0375
0376 particle.addPoint(layer1Pt);
0377
0378
0379 } else {
0380 reco::PFTrajectoryPoint dummy;
0381 particle.addPoint(dummy);
0382 }
0383
0384 if (fst.onLayer2()) {
0385 const RawParticle& rp = fst.layer2Entrance();
0386
0387 math::XYZPoint posLayer2(rp.x(), rp.y(), rp.z());
0388 math::XYZTLorentzVector momLayer2(rp.px(), rp.py(), rp.pz(), rp.e());
0389 reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2, posLayer2, momLayer2);
0390
0391 particle.addPoint(layer2Pt);
0392
0393
0394 } else {
0395 reco::PFTrajectoryPoint dummy;
0396 particle.addPoint(dummy);
0397 }
0398
0399 if (fst.onEcal()) {
0400 const RawParticle& rp = fst.ecalEntrance();
0401
0402 math::XYZPoint posECAL(rp.x(), rp.y(), rp.z());
0403 math::XYZTLorentzVector momECAL(rp.px(), rp.py(), rp.pz(), rp.e());
0404 reco::PFTrajectoryPoint ecalPt(-1, reco::PFTrajectoryPoint::ECALEntrance, posECAL, momECAL);
0405
0406 particle.addPoint(ecalPt);
0407
0408
0409 } else {
0410 reco::PFTrajectoryPoint dummy;
0411 particle.addPoint(dummy);
0412 }
0413
0414
0415 reco::PFTrajectoryPoint dummy;
0416 particle.addPoint(dummy);
0417
0418 if (fst.onHcal()) {
0419 const RawParticle& rpin = fst.hcalEntrance();
0420
0421 math::XYZPoint posHCALin(rpin.x(), rpin.y(), rpin.z());
0422 math::XYZTLorentzVector momHCALin(rpin.px(), rpin.py(), rpin.pz(), rpin.e());
0423 reco::PFTrajectoryPoint hcalPtin(-1, reco::PFTrajectoryPoint::HCALEntrance, posHCALin, momHCALin);
0424
0425 particle.addPoint(hcalPtin);
0426
0427 const RawParticle& rpout = fst.hcalExit();
0428
0429 math::XYZPoint posHCALout(rpout.x(), rpout.y(), rpout.z());
0430 math::XYZTLorentzVector momHCALout(rpout.px(), rpout.py(), rpout.pz(), rpout.e());
0431 reco::PFTrajectoryPoint hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit, posHCALout, momHCALout);
0432
0433 particle.addPoint(hcalPtout);
0434
0435 const RawParticle& rpho = fst.hoEntrance();
0436
0437 math::XYZPoint posHOEntrance(rpho.x(), rpho.y(), rpho.z());
0438 math::XYZTLorentzVector momHOEntrance(rpho.px(), rpho.py(), rpho.pz(), rpho.e());
0439 reco::PFTrajectoryPoint hoPtin(-1, reco::PFTrajectoryPoint::HOLayer, posHOEntrance, momHOEntrance);
0440
0441 particle.addPoint(hoPtin);
0442
0443 } else {
0444 reco::PFTrajectoryPoint dummy;
0445 particle.addPoint(dummy);
0446 }
0447
0448 pOutputPFSimParticleCollection->push_back(particle);
0449 }
0450
0451 iEvent.put(std::move(pOutputPFSimParticleCollection));
0452 }
0453
0454 LogDebug("PFSimParticleProducer") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
0455 }
0456
0457 void PFSimParticleProducer::getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID) {
0458 if (trackh.isValid()) {
0459 for (unsigned i = 0; i < trackh->size(); i++) {
0460 const reco::PFRecTrackRef ref(trackh, i);
0461
0462 for (auto const& hit : ref->trackRef()->recHits()) {
0463 if (hit->isValid()) {
0464 auto rechit = dynamic_cast<const FastTrackerRecHit*>(hit);
0465
0466 for (unsigned int st_index = 0; st_index < rechit->nSimTrackIds(); ++st_index) {
0467 recTrackSimID.push_back(rechit->simTrackId(st_index));
0468 }
0469 break;
0470 }
0471 }
0472 }
0473 }
0474 }