File indexing completed on 2023-03-17 11:26:45
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 #include <memory>
0025
0026
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034
0035 #include "DataFormats/Common/interface/Handle.h"
0036
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039
0040 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0041 #include "DataFormats/GeometrySurface/interface/Plane.h"
0042
0043 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0044 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0045
0046 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0047 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0048
0049 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0050
0051 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0052
0053 #include "MagneticField/Engine/interface/MagneticField.h"
0054 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0055
0056 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0057
0058 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0059 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0060 #include "SimDataFormats/Track/interface/SimTrack.h"
0061 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0062 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0063 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0064
0065 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0066 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0067 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0068
0069 #include "CLHEP/Random/RandFlat.h"
0070 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0071
0072
0073 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0074 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixStateInfo.h"
0075
0076 #include "TFile.h"
0077 #include "TTree.h"
0078
0079 #include <map>
0080
0081
0082
0083
0084
0085 class SteppingHelixPropagatorAnalyzer : public edm::one::EDAnalyzer<> {
0086 public:
0087 explicit SteppingHelixPropagatorAnalyzer(const edm::ParameterSet&);
0088 ~SteppingHelixPropagatorAnalyzer();
0089
0090 virtual void analyze(const edm::Event&, const edm::EventSetup&);
0091 virtual void endJob();
0092
0093 protected:
0094 struct GlobalSimHit {
0095 const PSimHit* hit;
0096 const Surface* surf;
0097 DetId id;
0098 CLHEP::Hep3Vector r3;
0099 CLHEP::Hep3Vector p3;
0100 };
0101
0102 void loadNtVars(int ind,
0103 int eType,
0104 int pStatus,
0105 int id,
0106 const CLHEP::Hep3Vector& p3,
0107 const CLHEP::Hep3Vector& r3,
0108 const CLHEP::Hep3Vector& p3R,
0109 const CLHEP::Hep3Vector& r3R,
0110 int charge,
0111 const AlgebraicSymMatrix66& cov);
0112
0113 FreeTrajectoryState getFromCLHEP(const CLHEP::Hep3Vector& p3,
0114 const CLHEP::Hep3Vector& r3,
0115 int charge,
0116 const AlgebraicSymMatrix66& cov,
0117 const MagneticField* field);
0118 void getFromFTS(const FreeTrajectoryState& fts,
0119 CLHEP::Hep3Vector& p3,
0120 CLHEP::Hep3Vector& r3,
0121 int& charge,
0122 AlgebraicSymMatrix66& cov);
0123
0124 void addPSimHits(const edm::Event& iEvent,
0125 const std::string instanceName,
0126 const edm::ESHandle<GlobalTrackingGeometry>& geom,
0127 std::vector<SteppingHelixPropagatorAnalyzer::GlobalSimHit>& hits) const;
0128
0129 private:
0130
0131 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
0132 const edm::ESGetToken<Propagator, TrackingComponentsRecord> propToken;
0133 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> globTrkGeomToken;
0134
0135 TFile* ntFile_;
0136 TTree* tr_;
0137
0138 int nPoints_;
0139 int q_[1000];
0140 int eType_[1000];
0141 int pStatus_[1000][3];
0142 float p3_[1000][9];
0143 float r3_[1000][9];
0144 int id_[1000];
0145 float p3R_[1000][3];
0146 float r3R_[1000][3];
0147 float covFlat_[1000][21];
0148
0149 bool debug_;
0150 int run_;
0151 int event_;
0152
0153 int trkIndOffset_;
0154
0155 bool doneMapping_;
0156
0157 bool radX0CorrectionMode_;
0158
0159 bool testPCAPropagation_;
0160
0161 bool ntupleTkHits_;
0162
0163 bool startFromPrevHit_;
0164
0165 std::string g4SimName_;
0166 edm::InputTag simTracksTag_;
0167 edm::InputTag simVertexesTag_;
0168 };
0169
0170
0171
0172
0173 SteppingHelixPropagatorAnalyzer::SteppingHelixPropagatorAnalyzer(const edm::ParameterSet& iConfig)
0174 : magFieldToken(esConsumes()),
0175 propToken(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0176 globTrkGeomToken(esConsumes()),
0177 simTracksTag_(iConfig.getParameter<edm::InputTag>("simTracksTag")),
0178 simVertexesTag_(iConfig.getParameter<edm::InputTag>("simVertexesTag")) {
0179
0180
0181 ntFile_ = new TFile(iConfig.getParameter<std::string>("NtFile").c_str(), "recreate");
0182 tr_ = new TTree("MuProp", "MuProp");
0183 tr_->Branch("nPoints", &nPoints_, "nPoints/I");
0184 tr_->Branch("q", q_, "q[nPoints]/I");
0185 tr_->Branch("pStatus", pStatus_, "pStatus[nPoints][3]/I");
0186 tr_->Branch("p3", p3_, "p3[nPoints][9]/F");
0187 tr_->Branch("r3", r3_, "r3[nPoints][9]/F");
0188 tr_->Branch("id", id_, "id[nPoints]/I");
0189 tr_->Branch("p3R", p3R_, "p3R[nPoints][3]/F");
0190 tr_->Branch("r3R", r3R_, "r3R[nPoints][3]/F");
0191 tr_->Branch("covFlat", covFlat_, "covFlat[nPoints][21]/F");
0192 tr_->Branch("run", &run_, "run/I");
0193 tr_->Branch("event_", &event_, "event/I");
0194
0195 trkIndOffset_ = iConfig.getParameter<int>("trkIndOffset");
0196 debug_ = iConfig.getParameter<bool>("debug");
0197 radX0CorrectionMode_ = iConfig.getParameter<bool>("radX0CorrectionMode");
0198
0199 testPCAPropagation_ = iConfig.getParameter<bool>("testPCAPropagation");
0200
0201 ntupleTkHits_ = iConfig.getParameter<bool>("ntupleTkHits");
0202
0203 startFromPrevHit_ = iConfig.getParameter<bool>("startFromPrevHit");
0204
0205 g4SimName_ = iConfig.getParameter<std::string>("g4SimName");
0206
0207
0208 for (int i = 0; i < 1000; ++i) {
0209 q_[i] = 0;
0210 eType_[i] = 0;
0211 for (int j = 0; j < 3; ++j)
0212 pStatus_[i][j] = 0;
0213 for (int j = 0; j < 9; ++j)
0214 p3_[i][j] = 0;
0215 for (int j = 0; j < 9; ++j)
0216 r3_[i][j] = 0;
0217 id_[i] = 0;
0218 for (int j = 0; j < 3; ++j)
0219 p3R_[i][j] = 0;
0220 for (int j = 0; j < 3; ++j)
0221 r3R_[i][j] = 0;
0222 for (int j = 0; j < 21; ++j)
0223 covFlat_[i][j] = 0;
0224 }
0225 }
0226
0227 SteppingHelixPropagatorAnalyzer::~SteppingHelixPropagatorAnalyzer() = default;
0228
0229
0230
0231
0232
0233
0234 void SteppingHelixPropagatorAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0235 constexpr char const* metname = "SteppingHelixPropagatorAnalyzer";
0236 (void)metname;
0237 using namespace edm;
0238 const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
0239 const ESHandle<Propagator> shProp = iSetup.getHandle(propToken);
0240 const ESHandle<Propagator> shPropAny = iSetup.getHandle(propToken);
0241 const ESHandle<GlobalTrackingGeometry> geomESH = iSetup.getHandle(globTrkGeomToken);
0242 if (debug_) {
0243 LogTrace(metname) << "Got GlobalTrackingGeometry " << std::endl;
0244 }
0245
0246 run_ = (int)iEvent.id().run();
0247 event_ = (int)iEvent.id().event();
0248 if (debug_) {
0249 LogTrace(metname) << "Begin for run:event ==" << run_ << ":" << event_ << std::endl;
0250 }
0251
0252 const double FPRP_MISMATCH = 150.;
0253 int pStatus = 0;
0254
0255 Handle<SimTrackContainer> simTracks;
0256 iEvent.getByLabel<SimTrackContainer>(simTracksTag_, simTracks);
0257 if (!simTracks.isValid()) {
0258 std::cout << "No tracks found" << std::endl;
0259 return;
0260 }
0261 if (debug_) {
0262 LogTrace(metname) << "Got simTracks of size " << simTracks->size() << std::endl;
0263 }
0264
0265 Handle<SimVertexContainer> simVertices;
0266 iEvent.getByLabel<SimVertexContainer>(simVertexesTag_, simVertices);
0267 if (!simVertices.isValid()) {
0268 std::cout << "No tracks found" << std::endl;
0269 return;
0270 }
0271 if (debug_) {
0272 LogTrace(metname) << "Got simVertices of size " << simVertices->size() << std::endl;
0273 }
0274
0275 std::vector<GlobalSimHit> allSimHits;
0276 allSimHits.clear();
0277
0278 addPSimHits(iEvent, "MuonDTHits", geomESH, allSimHits);
0279 addPSimHits(iEvent, "MuonCSCHits", geomESH, allSimHits);
0280 addPSimHits(iEvent, "MuonRPCHits", geomESH, allSimHits);
0281 addPSimHits(iEvent, "TrackerHitsPixelBarrelLowTof", geomESH, allSimHits);
0282 addPSimHits(iEvent, "TrackerHitsPixelEndcapLowTof", geomESH, allSimHits);
0283 addPSimHits(iEvent, "TrackerHitsTIBLowTof", geomESH, allSimHits);
0284 addPSimHits(iEvent, "TrackerHitsTIDLowTof", geomESH, allSimHits);
0285 addPSimHits(iEvent, "TrackerHitsTOBLowTof", geomESH, allSimHits);
0286 addPSimHits(iEvent, "TrackerHitsTECLowTof", geomESH, allSimHits);
0287
0288 SimTrackContainer::const_iterator tracksCI = simTracks->begin();
0289 for (; tracksCI != simTracks->end(); tracksCI++) {
0290 int trkPDG = tracksCI->type();
0291 if (abs(trkPDG) != 13) {
0292 if (debug_) {
0293 LogTrace(metname) << "Skip " << trkPDG << std::endl;
0294 }
0295 continue;
0296 }
0297 CLHEP::Hep3Vector p3T(tracksCI->momentum().x(), tracksCI->momentum().y(), tracksCI->momentum().z());
0298 if (p3T.mag() < 2.)
0299 continue;
0300
0301 int vtxInd = tracksCI->vertIndex();
0302 unsigned int trkInd = tracksCI->genpartIndex() - trkIndOffset_;
0303 CLHEP::Hep3Vector r3T(0., 0., 0.);
0304 if (vtxInd < 0) {
0305 std::cout << "Track with no vertex, defaulting to (0,0,0)" << std::endl;
0306 } else {
0307 r3T = CLHEP::Hep3Vector((*simVertices)[vtxInd].position().x(),
0308 (*simVertices)[vtxInd].position().y(),
0309 (*simVertices)[vtxInd].position().z());
0310 }
0311 AlgebraicSymMatrix66 covT = AlgebraicMatrixID();
0312 covT *= 1e-20;
0313
0314 CLHEP::Hep3Vector p3F, r3F;
0315 AlgebraicSymMatrix66 covF;
0316 int charge = trkPDG > 0 ? -1 : 1;
0317
0318 nPoints_ = 0;
0319 pStatus = 0;
0320 loadNtVars(nPoints_, 0, pStatus, 0, p3T, r3T, p3T, r3T, charge, covT);
0321 nPoints_++;
0322 FreeTrajectoryState ftsTrack = getFromCLHEP(p3T, r3T, charge, covT, &*bField);
0323 FreeTrajectoryState ftsStart = ftsTrack;
0324 SteppingHelixStateInfo siStart(ftsStart);
0325 SteppingHelixStateInfo siDest;
0326 TrajectoryStateOnSurface tSOSDest;
0327
0328 std::map<double, const GlobalSimHit*> simHitsByDistance;
0329 std::map<double, const GlobalSimHit*> simHitsByTof;
0330 for (std::vector<GlobalSimHit>::const_iterator allHitsCI = allSimHits.begin(); allHitsCI != allSimHits.end();
0331 allHitsCI++) {
0332 if (allHitsCI->hit->trackId() != trkInd)
0333 continue;
0334 if (abs(allHitsCI->hit->particleType()) != 13)
0335 continue;
0336 if (allHitsCI->p3.mag() < 0.5)
0337 continue;
0338
0339 double distance = (allHitsCI->r3 - r3T).mag();
0340 double tof = allHitsCI->hit->timeOfFlight();
0341 simHitsByDistance[distance] = &*allHitsCI;
0342 simHitsByTof[tof] = &*allHitsCI;
0343 }
0344
0345 {
0346 if (testPCAPropagation_) {
0347 FreeTrajectoryState ftsDest;
0348 GlobalPoint pDest1(10., 10., 0.);
0349 GlobalPoint pDest2(10., 10., 10.);
0350 const SteppingHelixPropagator* shPropAnyCPtr = dynamic_cast<const SteppingHelixPropagator*>(&*shPropAny);
0351
0352 ftsDest = shPropAnyCPtr->propagate(ftsStart, pDest1);
0353 std::cout << "----------------------------------------------" << std::endl;
0354 ftsDest = shPropAnyCPtr->propagate(ftsStart, pDest1, pDest2);
0355 std::cout << "----------------------------------------------" << std::endl;
0356 }
0357
0358
0359 std::map<double, const GlobalSimHit*>::const_iterator simHitsCI = simHitsByDistance.begin();
0360 for (; simHitsCI != simHitsByDistance.end(); simHitsCI++) {
0361 const GlobalSimHit* igHit = simHitsCI->second;
0362 const PSimHit* iHit = simHitsCI->second->hit;
0363
0364 if (debug_) {
0365 LogTrace(metname) << igHit->id.rawId() << " r3L:" << iHit->localPosition() << " r3G:" << igHit->r3
0366 << " p3L:" << iHit->momentumAtEntry() << " p3G:" << igHit->p3
0367 << " pId:" << iHit->particleType() << " tId:" << iHit->trackId() << std::endl;
0368 }
0369
0370 if (debug_) {
0371 LogTrace(metname) << "Will propagate to surface: " << igHit->surf->position() << " "
0372 << igHit->surf->rotation() << std::endl;
0373 }
0374 pStatus = 0;
0375 if (startFromPrevHit_) {
0376 if (simHitsCI != simHitsByDistance.begin()) {
0377 std::map<double, const GlobalSimHit*>::const_iterator simHitPrevCI = simHitsCI;
0378 simHitPrevCI--;
0379 const GlobalSimHit* gpHit = simHitPrevCI->second;
0380 ftsStart = getFromCLHEP(gpHit->p3, gpHit->r3, charge, covT, &*bField);
0381 siStart = SteppingHelixStateInfo(ftsTrack);
0382 }
0383 }
0384
0385 if (radX0CorrectionMode_) {
0386 const SteppingHelixPropagator* shPropCPtr = dynamic_cast<const SteppingHelixPropagator*>(&*shProp);
0387 shPropCPtr->propagate(siStart, *igHit->surf, siDest);
0388 if (siDest.isValid()) {
0389 siStart = siDest;
0390 siStart.getFreeState(ftsStart);
0391 getFromFTS(ftsStart, p3F, r3F, charge, covF);
0392 pStatus = 0;
0393 } else
0394 pStatus = 1;
0395 if (pStatus == 1 || (r3F - igHit->r3).mag() > FPRP_MISMATCH) {
0396
0397 siStart = SteppingHelixStateInfo(ftsTrack);
0398 pStatus = 1;
0399 }
0400 } else {
0401 tSOSDest = shProp->propagate(ftsStart, *igHit->surf);
0402 if (tSOSDest.isValid()) {
0403 ftsStart = *tSOSDest.freeState();
0404 getFromFTS(ftsStart, p3F, r3F, charge, covF);
0405 pStatus = 0;
0406 } else
0407 pStatus = 1;
0408 if (pStatus == 1 || (r3F - igHit->r3).mag() > FPRP_MISMATCH) {
0409
0410 ftsStart = ftsTrack;
0411 pStatus = 1;
0412 }
0413 }
0414
0415 if (debug_) {
0416 LogTrace(metname) << "Got to "
0417 << " r3Prp:" << r3F << " r3Hit:" << igHit->r3 << " p3Prp:" << p3F << " p3Hit:" << igHit->p3
0418 << " pPrp:" << p3F.mag() << " pHit:" << igHit->p3.mag() << std::endl;
0419 }
0420 loadNtVars(nPoints_, 0, pStatus, igHit->id.rawId(), p3F, r3F, igHit->p3, igHit->r3, charge, covF);
0421 nPoints_++;
0422 }
0423 }
0424 if (tr_)
0425 tr_->Fill();
0426 }
0427 }
0428
0429
0430
0431
0432 void SteppingHelixPropagatorAnalyzer::endJob() {
0433 ntFile_->cd();
0434 tr_->Write();
0435 delete ntFile_;
0436 ntFile_ = 0;
0437 }
0438
0439 void SteppingHelixPropagatorAnalyzer::loadNtVars(int ind,
0440 int eType,
0441 int pStatus,
0442 int id,
0443 const CLHEP::Hep3Vector& p3,
0444 const CLHEP::Hep3Vector& r3,
0445 const CLHEP::Hep3Vector& p3R,
0446 const CLHEP::Hep3Vector& r3R,
0447 int charge,
0448 const AlgebraicSymMatrix66& cov) {
0449 p3_[ind][eType * 3 + 0] = p3.x();
0450 p3_[ind][eType * 3 + 1] = p3.y();
0451 p3_[ind][eType * 3 + 2] = p3.z();
0452 r3_[ind][eType * 3 + 0] = r3.x();
0453 r3_[ind][eType * 3 + 1] = r3.y();
0454 r3_[ind][eType * 3 + 2] = r3.z();
0455 id_[ind] = id;
0456 p3R_[ind][0] = p3R.x();
0457 p3R_[ind][1] = p3R.y();
0458 p3R_[ind][2] = p3R.z();
0459 r3R_[ind][0] = r3R.x();
0460 r3R_[ind][1] = r3R.y();
0461 r3R_[ind][2] = r3R.z();
0462 int flatInd = 0;
0463 for (int i = 1; i <= cov.kRows; i++)
0464 for (int j = 1; j <= i; j++) {
0465 covFlat_[ind][flatInd] = cov(i - 1, j - 1);
0466 flatInd++;
0467 }
0468 q_[ind] = charge;
0469 eType_[ind] = eType;
0470 pStatus_[ind][eType] = pStatus;
0471 }
0472
0473 FreeTrajectoryState SteppingHelixPropagatorAnalyzer::getFromCLHEP(const CLHEP::Hep3Vector& p3,
0474 const CLHEP::Hep3Vector& r3,
0475 int charge,
0476 const AlgebraicSymMatrix66& cov,
0477 const MagneticField* field) {
0478 GlobalVector p3GV(p3.x(), p3.y(), p3.z());
0479 GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
0480 GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
0481
0482 CartesianTrajectoryError tCov(cov);
0483
0484 return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars);
0485 }
0486
0487 void SteppingHelixPropagatorAnalyzer::getFromFTS(const FreeTrajectoryState& fts,
0488 CLHEP::Hep3Vector& p3,
0489 CLHEP::Hep3Vector& r3,
0490 int& charge,
0491 AlgebraicSymMatrix66& cov) {
0492 GlobalVector p3GV = fts.momentum();
0493 GlobalPoint r3GP = fts.position();
0494
0495 p3.set(p3GV.x(), p3GV.y(), p3GV.z());
0496 r3.set(r3GP.x(), r3GP.y(), r3GP.z());
0497
0498 charge = fts.charge();
0499 cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
0500 }
0501
0502 void SteppingHelixPropagatorAnalyzer ::addPSimHits(
0503 const edm::Event& iEvent,
0504 const std::string instanceName,
0505 const edm::ESHandle<GlobalTrackingGeometry>& geom,
0506 std::vector<SteppingHelixPropagatorAnalyzer::GlobalSimHit>& hits) const {
0507 static std::string metname = "SteppingHelixPropagatorAnalyzer";
0508 edm::Handle<edm::PSimHitContainer> handle;
0509 iEvent.getByLabel(g4SimName_, instanceName, handle);
0510 if (!handle.isValid()) {
0511 std::cout << "No hits found" << std::endl;
0512 return;
0513 }
0514 if (debug_) {
0515 LogTrace(metname) << "Got " << instanceName << " of size " << handle->size() << std::endl;
0516 }
0517
0518 edm::PSimHitContainer::const_iterator pHits_CI = handle->begin();
0519 for (; pHits_CI != handle->end(); pHits_CI++) {
0520 int dtId = pHits_CI->detUnitId();
0521 DetId wId(dtId);
0522
0523 const GeomDet* layer = geom->idToDet(wId);
0524
0525 GlobalSimHit gHit;
0526 gHit.hit = &*pHits_CI;
0527 gHit.surf = &layer->surface();
0528 gHit.id = wId;
0529
0530
0531 float dzFrac = gHit.hit->entryPoint().z() / (gHit.hit->exitPoint().z() - gHit.hit->entryPoint().z());
0532 GlobalPoint r3Hit =
0533 gHit.surf->toGlobal(gHit.hit->entryPoint() - dzFrac * (gHit.hit->exitPoint() - gHit.hit->entryPoint()));
0534 gHit.r3.set(r3Hit.x(), r3Hit.y(), r3Hit.z());
0535 GlobalVector p3Hit = gHit.surf->toGlobal(gHit.hit->momentumAtEntry());
0536 gHit.p3.set(p3Hit.x(), p3Hit.y(), p3Hit.z());
0537 hits.push_back(gHit);
0538 }
0539 }
0540
0541
0542 DEFINE_FWK_MODULE(SteppingHelixPropagatorAnalyzer);