File indexing completed on 2022-01-25 05:52:16
0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h" //For define_fwk_module
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010
0011
0012
0013
0014
0015 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0016 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0017 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0018 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0021 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0022 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0023
0024
0025 #include "MagneticField/Engine/interface/MagneticField.h"
0026 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0027 #include "SimG4Core/MagneticField/interface/Field.h"
0028 #include "SimG4Core/MagneticField/interface/FieldBuilder.h"
0029
0030
0031 #include "TrackPropagation/Geant4e/interface/ConvertFromToCLHEP.h"
0032 #include "TrackPropagation/Geant4e/interface/Geant4ePropagator.h"
0033 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0034 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0035
0036
0037 #include "SimDataFormats/Track/interface/SimTrack.h"
0038 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0039 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0040 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0041 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0042 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0043
0044
0045 #include "G4TransportationManager.hh"
0046
0047
0048 #include "TFile.h"
0049 #include "TH1.h"
0050 #include "TH2.h"
0051
0052 using namespace std;
0053
0054 enum testMuChamberType { DT, RPC, CSC };
0055
0056 class Geant4ePropagatorAnalyzer : public edm::one::EDAnalyzer<> {
0057 public:
0058 explicit Geant4ePropagatorAnalyzer(const edm::ParameterSet &);
0059 ~Geant4ePropagatorAnalyzer() override {}
0060
0061 void analyze(const edm::Event &, const edm::EventSetup &) override;
0062 void endJob() override;
0063 void beginJob() override;
0064 void iterateOverHits(edm::Handle<edm::PSimHitContainer> simHits,
0065 testMuChamberType muonChamberType,
0066 unsigned int trkIndex,
0067 const FreeTrajectoryState &ftsTrack);
0068
0069 protected:
0070
0071 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
0072 const edm::ESGetToken<DTGeometry, MuonGeometryRecord> dtGeomToken;
0073 const edm::ESGetToken<RPCGeometry, MuonGeometryRecord> rpcGeomToken;
0074 const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> cscGeomToken;
0075
0076 int theRun;
0077 int theEvent;
0078
0079 Propagator *thePropagator;
0080
0081
0082
0083 edm::ESHandle<DTGeometry> theDTGeomESH;
0084 edm::ESHandle<RPCGeometry> theRPCGeomESH;
0085 edm::ESHandle<CSCGeometry> theCSCGeomESH;
0086
0087
0088 int fStudyStation;
0089
0090
0091 float fBeamCenter;
0092 float fBeamInterval;
0093
0094
0095 TFile *theRootFile;
0096
0097 TH1F *fDistanceSHLayer;
0098
0099 TH1F *fDistance;
0100 TH1F *fDistanceSt[4];
0101
0102 TH1F *fHitR;
0103 TH1F *fHitRho;
0104 TH1F *fHitEta;
0105 TH1F *fHitPhi;
0106 TH1F *fHitPhi1L;
0107 TH2F *fHitRVsPhi;
0108
0109 TH1F *fLayerR;
0110 TH1F *fLayerRho;
0111 TH1F *fLayerEta;
0112 TH1F *fLayerPhi;
0113 TH2F *fLayerRVsPhi;
0114
0115 TH1F *fExtrapR;
0116 TH1F *fExtrapRho;
0117 TH1F *fExtrapEta;
0118 TH1F *fExtrapPhi;
0119 TH2F *fExtrapRVsPhi;
0120
0121 TH1F *fDeltaRo;
0122 TH1F *fDeltaEta;
0123 TH1F *fDeltaPhi;
0124
0125
0126 TH1F *fStationPosPhi;
0127 TH1F *fSectorPosPhi;
0128 TH1F *fSLayerPosPhi;
0129 TH1F *fLayerPosPhi;
0130
0131 TH1F *fStationNegPhi;
0132 TH1F *fSectorNegPhi;
0133 TH1F *fSLayerNegPhi;
0134 TH1F *fLayerNegPhi;
0135
0136 edm::InputTag G4VtxSrc_;
0137 edm::InputTag G4TrkSrc_;
0138 };
0139
0140 Geant4ePropagatorAnalyzer::Geant4ePropagatorAnalyzer(const edm::ParameterSet &iConfig)
0141 : magFieldToken(esConsumes()),
0142 dtGeomToken(esConsumes()),
0143 rpcGeomToken(esConsumes()),
0144 cscGeomToken(esConsumes()),
0145 theRun(-1),
0146 theEvent(-1),
0147 thePropagator(nullptr),
0148 G4VtxSrc_(iConfig.getParameter<edm::InputTag>("G4VtxSrc")),
0149 G4TrkSrc_(iConfig.getParameter<edm::InputTag>("G4TrkSrc")) {
0150
0151 fStudyStation = iConfig.getParameter<int>("StudyStation");
0152
0153 fBeamCenter = iConfig.getParameter<double>("BeamCenter");
0154 fBeamInterval = iConfig.getParameter<double>("BeamInterval");
0155
0156
0157
0158
0159
0160
0161 theRootFile = new TFile(iConfig.getParameter<std::string>("RootFile").c_str(), "recreate");
0162
0163
0164 fDistanceSHLayer = new TH1F("fDistanceSHLayer", "Distance(sim hit - layer)", 200, -1, 1);
0165
0166
0167 fDistance = new TH1F("fDistance", "Distance(sim hit - extrap)", 150, 0, 300);
0168
0169 fDistanceSt[0] = new TH1F("fDistance_St1", "Distance(sim hit - extrap) for station 1", 150, 0, 300);
0170 fDistanceSt[1] = new TH1F("fDistance_St2", "Distance(sim hit - extrap) for station 2", 150, 0, 300);
0171 fDistanceSt[2] = new TH1F("fDistance_St3", "Distance(sim hit - extrap) for station 3", 150, 0, 300);
0172 fDistanceSt[3] = new TH1F("fDistance_St4", "Distance(sim hit - extrap) for station 4", 150, 0, 300);
0173
0174
0175 fHitR = new TH1F("fHitR", "R^{sim hit}", 300, 400, 1000);
0176 fHitRho = new TH1F("fHitRho", "#rho^{sim hit}", 300, 400, 1000);
0177 fHitEta = new TH1F("fHitEta", "#eta^{sim hit}", 100, -0.1, 0.1);
0178 fHitPhi = new TH1F("fHitPhi", "#varphi^{sim hit}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0179 fHitPhi1L = new TH1F(
0180 "fHitPhi1L", "#varphi^{sim hit} for 1st layer", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0181
0182 fHitRVsPhi = new TH2F("fHitRVsPhi",
0183 "R vs. #varphi for sim hits",
0184 60,
0185 400,
0186 1000,
0187 80,
0188 fBeamCenter - fBeamInterval,
0189 fBeamCenter + fBeamInterval);
0190
0191
0192 fLayerR = new TH1F("fLayerR", "R^{layer}", 300, 400, 1000);
0193 fLayerRho = new TH1F("fLayerRho", "#rho^{layer}", 300, 400, 1000);
0194 fLayerEta = new TH1F("fLayerEta", "#eta^{layer}", 100, -0.1, 0.1);
0195 fLayerPhi = new TH1F("fLayerPhi", "#varphi^{layer}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0196 fLayerRVsPhi = new TH2F("fLayerRVsPhi",
0197 "R vs. #varphi for layers",
0198 60,
0199 400,
0200 1000,
0201 80,
0202 fBeamCenter - fBeamInterval,
0203 fBeamCenter + fBeamInterval);
0204
0205
0206 fExtrapR = new TH1F("fExtrapR", "R^{extrap. hit}", 300, 400, 1000);
0207 fExtrapRho = new TH1F("fExtrapRho", "#rho^{extrap. hit}", 300, 400, 1000);
0208 fExtrapEta = new TH1F("fExtrapEta", "#eta^{extrap. hit}", 100, -0.1, 0.1);
0209 fExtrapPhi =
0210 new TH1F("fExtrapPhi", "#varphi^{extrap. hit}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0211
0212 fExtrapRVsPhi = new TH2F("fExtrapRVsPhi",
0213 "R vs. #varphi for extrap. hits",
0214 60,
0215 400,
0216 1000,
0217 80,
0218 fBeamCenter - fBeamInterval,
0219 fBeamCenter + fBeamInterval);
0220
0221
0222 fDeltaRo = new TH1F("fDeltaRo", "#Delta(#rho^{sim}, #rho^{extrap})", 100, 0, 200);
0223 fDeltaEta = new TH1F("fDeltaEta", "#Delta(#eta^{sim}, #eta^{extrap})", 100, -1, 1);
0224 fDeltaPhi = new TH1F("fDeltaPhi", "#Delta(#varphi^{sim}, #varphi^{extrap})", 120, -30, 30);
0225
0226
0227 fStationPosPhi = new TH1F("fStationPosPhi", "Station with positive #varphi", 6, -0.5, 5.5);
0228 fSectorPosPhi = new TH1F("fSectorPosPhi", "Sector with positive #varphi", 6, -0.5, 5.5);
0229 fSLayerPosPhi = new TH1F("fSLayerPosPhi", "Superlayer with positive #varphi", 6, -0.5, 5.5);
0230 fLayerPosPhi = new TH1F("fLayerPosPhi", "Layer with positive #varphi", 6, -0.5, 5.5);
0231 fStationNegPhi = new TH1F("fStationNegPhi", "Station with negative #varphi", 6, -0.5, 5.5);
0232 fSectorNegPhi = new TH1F("fSectorNegPhi", "Sector with negative #varphi", 6, -0.5, 5.5);
0233 fSLayerNegPhi = new TH1F("fSLayerNegPhi", "Superlayer with negative #varphi", 6, -0.5, 5.5);
0234 fLayerNegPhi = new TH1F("fLayerNegPhi", "Layer with negative #varphi", 6, -0.5, 5.5);
0235
0236
0237
0238 }
0239
0240 void Geant4ePropagatorAnalyzer::beginJob() { LogDebug("Geant4e") << "Nothing done in beginJob..."; }
0241
0242 void Geant4ePropagatorAnalyzer::endJob() {
0243 fDistanceSHLayer->Write();
0244 fDistance->Write();
0245 for (unsigned int i = 0; i < 4; i++)
0246 fDistanceSt[i]->Write();
0247
0248 fHitR->Write();
0249 fHitRho->Write();
0250 fHitEta->Write();
0251 fHitPhi->Write();
0252 fHitPhi1L->Write();
0253 fHitRVsPhi->Write();
0254
0255 fLayerR->Write();
0256 fLayerRho->Write();
0257 fLayerEta->Write();
0258 fLayerPhi->Write();
0259 fLayerRVsPhi->Write();
0260
0261 fExtrapR->Write();
0262 fExtrapRho->Write();
0263 fExtrapEta->Write();
0264 fExtrapPhi->Write();
0265 fExtrapRVsPhi->Write();
0266
0267 fDeltaRo->Write();
0268 fDeltaEta->Write();
0269 fDeltaPhi->Write();
0270
0271
0272 fStationPosPhi->Write();
0273 fSectorPosPhi->Write();
0274 fSLayerPosPhi->Write();
0275 fLayerPosPhi->Write();
0276 fStationNegPhi->Write();
0277 fSectorNegPhi->Write();
0278 fSLayerNegPhi->Write();
0279 fLayerNegPhi->Write();
0280
0281 theRootFile->Close();
0282
0283 }
0284
0285 void Geant4ePropagatorAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0286 using namespace edm;
0287
0288 LogDebug("Geant4e") << "Starting analyze...";
0289
0290
0291
0292 const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
0293 if (bField.isValid())
0294 LogDebug("Geant4e") << "G4e -- Magnetic field is valid. Value in (0,0,0): "
0295 << bField->inTesla(GlobalPoint(0, 0, 0)).mag() << " Tesla";
0296 else
0297 LogError("Geant4e") << "G4e -- NO valid Magnetic field";
0298
0299
0300
0301
0302
0303 theDTGeomESH = iSetup.getHandle(dtGeomToken);
0304 LogDebug("Geant4e") << "Got DTGeometry";
0305
0306
0307 theCSCGeomESH = iSetup.getHandle(cscGeomToken);
0308 LogDebug("Geant4e") << "Got CSCGeometry";
0309
0310
0311 theRPCGeomESH = iSetup.getHandle(rpcGeomToken);
0312 LogDebug("Geant4e") << "Got RPCGeometry";
0313
0314
0315
0316 theRun = (int)iEvent.id().run();
0317 theEvent = (int)iEvent.id().event();
0318 LogDebug("Geant4e") << "G4e -- Begin for run/event ==" << theRun << "/" << theEvent
0319 << " ---------------------------------";
0320
0321
0322
0323 if (!thePropagator)
0324 thePropagator = new Geant4ePropagator(&*bField);
0325
0326 if (thePropagator)
0327 LogDebug("Geant4e") << "Propagator built!";
0328 else
0329 LogError("Geant4e") << "Could not build propagator!";
0330
0331
0332
0333 Handle<SimTrackContainer> simTracks;
0334 iEvent.getByLabel<SimTrackContainer>(G4TrkSrc_, simTracks);
0335 if (!simTracks.isValid()) {
0336 LogWarning("Geant4e") << "No tracks found" << std::endl;
0337 return;
0338 }
0339 LogDebug("Geant4e") << "G4e -- Got simTracks of size " << simTracks->size();
0340
0341 Handle<SimVertexContainer> simVertices;
0342 iEvent.getByLabel<SimVertexContainer>(G4VtxSrc_, simVertices);
0343 if (!simVertices.isValid()) {
0344 LogWarning("Geant4e") << "No tracks found" << std::endl;
0345 return;
0346 }
0347 LogDebug("Geant4e") << "Got simVertices of size " << simVertices->size();
0348
0349
0350
0351 Handle<PSimHitContainer> simHitsDT;
0352 iEvent.getByLabel("g4SimHits", "MuonDTHits", simHitsDT);
0353 if (!simHitsDT.isValid()) {
0354 LogWarning("Geant4e") << "No hits found" << std::endl;
0355 return;
0356 }
0357 LogDebug("Geant4e") << "Got MuonDTHits of size " << simHitsDT->size();
0358
0359 Handle<PSimHitContainer> simHitsCSC;
0360 iEvent.getByLabel("g4SimHits", "MuonCSCHits", simHitsCSC);
0361 if (!simHitsCSC.isValid()) {
0362 LogWarning("Geant4e") << "No hits found" << std::endl;
0363 return;
0364 }
0365 LogDebug("Geant4e") << "Got MuonCSCHits of size " << simHitsCSC->size();
0366
0367 Handle<PSimHitContainer> simHitsRPC;
0368 iEvent.getByLabel("g4SimHits", "MuonRPCHits", simHitsRPC);
0369 if (!simHitsRPC.isValid()) {
0370 LogWarning("Geant4e") << "No hits found" << std::endl;
0371 return;
0372 }
0373 LogDebug("Geant4e") << "Got MuonRPCHits of size " << simHitsRPC->size();
0374
0375
0376
0377
0378
0379 unsigned int counter = 0;
0380
0381 for (SimTrackContainer::const_iterator simTracksIt = simTracks->begin(); simTracksIt != simTracks->end();
0382 simTracksIt++) {
0383
0384 counter++;
0385 LogDebug("Geant4e") << "G4e -- Iterating over " << counter << "rd track. Number: " << simTracksIt->genpartIndex()
0386 << "------------------";
0387
0388
0389 int simTrackPDG = simTracksIt->type();
0390 if (abs(simTrackPDG) != 13) {
0391 continue;
0392 }
0393 LogDebug("Geant4e") << "G4e -- Track PDG " << simTrackPDG;
0394
0395
0396
0397
0398
0399 int trkPDG = simTracksIt->type();
0400 if (abs(trkPDG) != 13) {
0401 LogDebug("Geant4e") << "Track is not a muon: " << trkPDG;
0402 continue;
0403 } else
0404 LogDebug("Geant4e") << "Found a muon track " << trkPDG;
0405
0406
0407 GlobalVector p3T = TrackPropagation::hep3VectorToGlobalVector(
0408 CLHEP::Hep3Vector(simTracksIt->momentum().x(), simTracksIt->momentum().y(), simTracksIt->momentum().z()));
0409 if (p3T.perp() < 2.) {
0410 LogDebug("Geant4e") << "Track PT is too low: " << p3T.perp();
0411 continue;
0412 } else {
0413 LogDebug("Geant4e") << "*** Phi (rad): " << p3T.phi() << " - Phi(deg)" << p3T.phi().degrees();
0414 LogDebug("Geant4e") << "Track PT is enough.";
0415 LogDebug("Geant4e") << "Track P.: " << p3T << "\nTrack P.: PT=" << p3T.perp() << "\tEta=" << p3T.eta()
0416 << "\tPhi=" << p3T.phi().degrees() << "--> Rad: Phi=" << p3T.phi();
0417 }
0418
0419
0420 int vtxInd = simTracksIt->vertIndex();
0421 GlobalPoint r3T(0., 0., 0.);
0422 if (vtxInd < 0)
0423 LogDebug("Geant4e") << "Track with no vertex, defaulting to (0,0,0)";
0424 else
0425
0426 r3T = TrackPropagation::hep3VectorToGlobalPoint(CLHEP::Hep3Vector((*simVertices)[vtxInd].position().x(),
0427 (*simVertices)[vtxInd].position().y(),
0428 (*simVertices)[vtxInd].position().z()));
0429
0430 LogDebug("Geant4e") << "Init point: " << r3T << "\nInit point Ro=" << r3T.perp() << "\tEta=" << r3T.eta()
0431 << "\tPhi=" << r3T.phi().degrees();
0432
0433
0434 int charge = trkPDG > 0 ? -1 : 1;
0435 LogDebug("Geant4e") << "Track charge = " << charge;
0436
0437
0438 CurvilinearTrajectoryError covT;
0439
0440
0441
0442 GlobalTrajectoryParameters trackPars(r3T, p3T, charge, &*bField);
0443 FreeTrajectoryState ftsTrack(trackPars, covT);
0444
0445
0446 unsigned int trkInd = simTracksIt->genpartIndex();
0447
0448
0449
0450 iterateOverHits(simHitsDT, DT, trkInd, ftsTrack);
0451
0452
0453
0454
0455
0456
0457
0458 }
0459 }
0460
0461 void Geant4ePropagatorAnalyzer::iterateOverHits(edm::Handle<edm::PSimHitContainer> simHits,
0462 testMuChamberType muonChamberType,
0463 unsigned int trkIndex,
0464 const FreeTrajectoryState &ftsTrack) {
0465 using namespace edm;
0466
0467 if (muonChamberType == DT)
0468 LogDebug("Geant4e") << "G4e -- Iterating over DT hits...";
0469 else if (muonChamberType == RPC)
0470 LogDebug("Geant4e") << "G4e -- Iterating over RPC hits...";
0471 else if (muonChamberType == CSC)
0472 LogDebug("Geant4e") << "G4e -- Iterating over CSC hits...";
0473
0474 for (PSimHitContainer::const_iterator simHitIt = simHits->begin(); simHitIt != simHits->end(); simHitIt++) {
0475
0476
0477 if (simHitIt->trackId() != trkIndex) {
0478 LogDebug("Geant4e") << "Hit (in tr " << simHitIt->trackId() << ") does not belong to track " << trkIndex;
0479 continue;
0480 }
0481
0482 LogDebug("Geant4e") << "G4e -- Hit belongs to track " << trkIndex;
0483
0484
0485
0486 int trkPDG = simHitIt->particleType();
0487 if (abs(trkPDG) != 13) {
0488 LogDebug("Geant4e") << "Associated track is not a muon: " << trkPDG;
0489 continue;
0490 }
0491 LogDebug("Geant4e") << "G4e -- Found a hit corresponding to a muon " << trkPDG;
0492
0493
0494
0495
0496 const GeomDet *layer = nullptr;
0497
0498 DTWireId *wIdDT = nullptr;
0499 if (muonChamberType == DT) {
0500 wIdDT = new DTWireId(simHitIt->detUnitId());
0501 int station = wIdDT->station();
0502 LogDebug("Geant4e") << "G4e -- DT Chamber. Station: " << station << ". DetId: " << wIdDT->layerId();
0503
0504 if (fStudyStation != -1 && fStudyStation != station)
0505 continue;
0506
0507 layer = theDTGeomESH->layer(*wIdDT);
0508 if (layer == nullptr) {
0509 LogDebug("Geant4e") << "Failed to get detector unit";
0510 continue;
0511 }
0512 }
0513
0514 else if (muonChamberType == RPC) {
0515 RPCDetId rpcId(simHitIt->detUnitId());
0516 layer = theRPCGeomESH->idToDet(rpcId);
0517 if (layer == nullptr) {
0518 LogDebug("Geant4e") << "Failed to get detector unit";
0519 continue;
0520 }
0521 }
0522
0523
0524 else if (muonChamberType == CSC) {
0525 CSCDetId cscId(simHitIt->detUnitId());
0526 layer = theCSCGeomESH->idToDet(cscId);
0527 if (layer == nullptr) {
0528 LogDebug("Geant4e") << "Failed to get detector unit";
0529 continue;
0530 }
0531 }
0532
0533 const Surface &surf = layer->surface();
0534 if (layer->geographicalId() != DTLayerId(simHitIt->detUnitId()))
0535 LogError("Geant4e") << "ERROR: wrong DetId";
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547 GlobalVector p3Hit = surf.toGlobal(simHitIt->momentumAtEntry());
0548 if (p3Hit.perp() < 0.5)
0549 continue;
0550 GlobalPoint posHit = surf.toGlobal(simHitIt->localPosition());
0551 Point3DBase<float, GlobalTag> surfpos = surf.position();
0552 LogDebug("Geant4e") << "G4e -- Layer position: " << surfpos << " cm"
0553 << "\nG4e -- Ro=" << surfpos.perp() << "\tEta=" << surfpos.eta()
0554 << "\tPhi=" << surfpos.phi().degrees() << "deg";
0555 LogDebug("Geant4e") << "G4e -- Sim Hit position: " << posHit << " cm"
0556 << "\nG4e -- Ro=" << posHit.perp() << "\tEta=" << posHit.eta()
0557 << "\tPhi=" << posHit.phi().degrees() << "deg"
0558 << "\t localpos=" << simHitIt->localPosition();
0559
0560 const Plane *bp = dynamic_cast<const Plane *>(&surf);
0561 if (bp != nullptr) {
0562 Float_t distance = bp->localZ(posHit);
0563 LogDebug("Geant4e") << "\nG4e -- Distance from plane to sim hit: " << distance << "cm";
0564 fDistanceSHLayer->Fill(distance);
0565 } else {
0566 LogWarning("Geant4e") << "G4e -- Layer is not a Plane!!!";
0567 fDistanceSHLayer->Fill(-1);
0568 }
0569
0570 LogDebug("Geant4e") << "Sim Hit Momentum PT=" << p3Hit.perp() << "\tEta=" << p3Hit.eta()
0571 << "\tPhi=" << p3Hit.phi().degrees() << "deg";
0572
0573
0574
0575 TrajectoryStateOnSurface tSOSDest = thePropagator->propagate(ftsTrack, surf);
0576
0577
0578
0579 GlobalPoint posExtrap = tSOSDest.freeState()->position();
0580
0581
0582
0583
0584 GlobalVector posDistance = posExtrap - posHit;
0585 float distance = posDistance.mag();
0586 float simhitphi = posHit.phi().degrees();
0587 float layerphi = surfpos.phi().degrees();
0588 float extrapphi = posExtrap.phi().degrees();
0589 LogDebug("Geant4e") << "G4e -- Difference between hit and final position: " << distance << " cm\n"
0590 << "G4e -- Transversal difference between hit and final position: " << posDistance.perp()
0591 << " cm.";
0592 LogDebug("Geant4e") << "G4e -- Extrapolated position:" << posExtrap << " cm\n"
0593 << "G4e -- (Rho, eta, phi): (" << posExtrap.perp() << " cm, " << posExtrap.eta() << ", "
0594 << extrapphi << " deg)";
0595 LogDebug("Geant4e") << "G4e -- Hit position: " << posHit << " cm\n"
0596 << "G4e -- (Rho, eta, phi): (" << posHit.perp() << " cm, " << posHit.eta() << ", "
0597 << simhitphi << " deg)";
0598
0599 fDistance->Fill(distance);
0600 fDistanceSt[wIdDT->station() - 1]->Fill(distance);
0601
0602 fHitR->Fill(posHit.mag());
0603 fHitRho->Fill(posHit.perp());
0604 fHitEta->Fill(posHit.eta());
0605 fHitPhi->Fill(simhitphi);
0606 if (posHit.perp() < 500)
0607 fHitPhi1L->Fill(simhitphi);
0608 fHitRVsPhi->Fill(posHit.mag(), simhitphi);
0609
0610 fLayerR->Fill(surfpos.mag());
0611 fLayerRho->Fill(surfpos.perp());
0612 fLayerEta->Fill(surfpos.eta());
0613 fLayerPhi->Fill(layerphi);
0614 fLayerRVsPhi->Fill(surfpos.mag(), layerphi);
0615
0616 fExtrapR->Fill(posExtrap.mag());
0617 fExtrapRho->Fill(posExtrap.perp());
0618 fExtrapEta->Fill(posExtrap.eta());
0619 fExtrapPhi->Fill(extrapphi);
0620 fExtrapRVsPhi->Fill(posExtrap.mag(), extrapphi);
0621
0622 fDeltaRo->Fill(posExtrap.perp() - posHit.perp());
0623 fDeltaEta->Fill(posExtrap.eta() - posHit.eta());
0624 fDeltaPhi->Fill(extrapphi - simhitphi);
0625
0626 if (wIdDT) {
0627 if (simhitphi > 0) {
0628 fStationPosPhi->Fill(wIdDT->station());
0629 fSectorPosPhi->Fill(wIdDT->sector());
0630 fSLayerPosPhi->Fill(wIdDT->superlayer());
0631 fLayerPosPhi->Fill(wIdDT->layer());
0632 } else {
0633 fStationNegPhi->Fill(wIdDT->station());
0634 fSectorNegPhi->Fill(wIdDT->sector());
0635 fSLayerNegPhi->Fill(wIdDT->superlayer());
0636 fLayerNegPhi->Fill(wIdDT->layer());
0637 }
0638
0639 delete wIdDT;
0640 }
0641
0642 }
0643 }
0644
0645
0646 DEFINE_FWK_MODULE(Geant4ePropagatorAnalyzer);