Geant4e/test/Geant4ePropagatorAnalyzer.cc

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 //- Timing
0012 //#include "Utilities/Timing/interface/TimingReport.h"
0013
0014 //- Geometry
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 //- Magnetic field
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 //- Propagator
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 //- SimHits, Tracks and Vertices
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 //- Geant4
0045 #include "G4TransportationManager.hh"
0046
0047 //- ROOT
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 = default;
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   // event setup tokens
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   // std::auto_ptr<sim::FieldBuilder> theFieldBuilder;
0081
0082   // Geometry
0083   edm::ESHandle<DTGeometry> theDTGeomESH;    // DTs
0084   edm::ESHandle<RPCGeometry> theRPCGeomESH;  // RPC
0085   edm::ESHandle<CSCGeometry> theCSCGeomESH;  // CSC
0086
0087   // Station that we want to study. -1 for all
0088   int fStudyStation;
0089
0090   // Single muons Beam direction (phi) and interval to plot simhits
0091   float fBeamCenter;
0092   float fBeamInterval;
0093
0094   // Histograms and ROOT stuff
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   // Studies on Phi distribution
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   // event data tokens
0137   edm::InputTag G4TrkSrc_;
0138   edm::InputTag G4VtxSrc_;
0139   const edm::EDGetTokenT<edm::SimTrackContainer> simTrackToken_;
0140   const edm::EDGetTokenT<edm::SimVertexContainer> simVertexToken_;
0141   const edm::EDGetTokenT<edm::PSimHitContainer> simHitsDTToken_;
0142   const edm::EDGetTokenT<edm::PSimHitContainer> simHitsCSCToken_;
0143   const edm::EDGetTokenT<edm::PSimHitContainer> simHitsRPCToken_;
0144 };
0145
0146 Geant4ePropagatorAnalyzer::Geant4ePropagatorAnalyzer(const edm::ParameterSet &iConfig)
0147     : magFieldToken(esConsumes()),
0148       dtGeomToken(esConsumes()),
0149       rpcGeomToken(esConsumes()),
0150       cscGeomToken(esConsumes()),
0151       theRun(-1),
0152       theEvent(-1),
0153       thePropagator(nullptr),
0154       G4TrkSrc_(iConfig.getParameter<edm::InputTag>("G4TrkSrc")),
0155       G4VtxSrc_(iConfig.getParameter<edm::InputTag>("G4VtxSrc")),
0156       simTrackToken_(consumes<edm::SimTrackContainer>(G4TrkSrc_)),
0157       simVertexToken_(consumes<edm::SimVertexContainer>(G4VtxSrc_)),
0158       simHitsDTToken_(consumes<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "MuonDTHits"))),
0159       simHitsCSCToken_(consumes<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "MuonCSCHits"))),
0160       simHitsRPCToken_(consumes<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "MuonRPCHits"))) {
0161   // debug_ = iConfig.getParameter<bool>("debug");
0162   fStudyStation = iConfig.getParameter<int>("StudyStation");
0163
0164   fBeamCenter = iConfig.getParameter<double>("BeamCenter");
0165   fBeamInterval = iConfig.getParameter<double>("BeamInterval");
0166
0167   ///////////////////////////////////////////////////////////////////////////////////
0168   // Histograms
0169   //
0170
0171   // Output ROOT file
0172   theRootFile = new TFile(iConfig.getParameter<std::string>("RootFile").c_str(), "recreate");
0173
0174   // Distance between Sim Hit and associated Layer
0175   fDistanceSHLayer = new TH1F("fDistanceSHLayer", "Distance(sim hit - layer)", 200, -1, 1);
0176
0177   // Distance between simhit and extrapolation
0178   fDistance = new TH1F("fDistance", "Distance(sim hit - extrap)", 150, 0, 300);
0179   // Distance between simhit and extrapolation
0180   fDistanceSt[0] = new TH1F("fDistance_St1", "Distance(sim hit - extrap) for station 1", 150, 0, 300);
0181   fDistanceSt[1] = new TH1F("fDistance_St2", "Distance(sim hit - extrap) for station 2", 150, 0, 300);
0182   fDistanceSt[2] = new TH1F("fDistance_St3", "Distance(sim hit - extrap) for station 3", 150, 0, 300);
0183   fDistanceSt[3] = new TH1F("fDistance_St4", "Distance(sim hit - extrap) for station 4", 150, 0, 300);
0184
0185   // Simulated hits
0186   fHitR = new TH1F("fHitR", "R^{sim hit}", 300, 400, 1000);
0187   fHitRho = new TH1F("fHitRho", "#rho^{sim hit}", 300, 400, 1000);
0188   fHitEta = new TH1F("fHitEta", "#eta^{sim hit}", 100, -0.1, 0.1);
0189   fHitPhi = new TH1F("fHitPhi", "#varphi^{sim hit}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0190   fHitPhi1L = new TH1F(
0191       "fHitPhi1L", "#varphi^{sim hit} for 1st layer", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0192
0193   fHitRVsPhi = new TH2F("fHitRVsPhi",
0194                         "R vs. #varphi for sim hits",
0195                         60,
0196                         400,
0197                         1000,
0198                         80,
0199                         fBeamCenter - fBeamInterval,
0200                         fBeamCenter + fBeamInterval);
0201
0202   // Position of layers
0203   fLayerR = new TH1F("fLayerR", "R^{layer}", 300, 400, 1000);
0204   fLayerRho = new TH1F("fLayerRho", "#rho^{layer}", 300, 400, 1000);
0205   fLayerEta = new TH1F("fLayerEta", "#eta^{layer}", 100, -0.1, 0.1);
0206   fLayerPhi = new TH1F("fLayerPhi", "#varphi^{layer}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0207   fLayerRVsPhi = new TH2F("fLayerRVsPhi",
0208                           "R vs. #varphi for layers",
0209                           60,
0210                           400,
0211                           1000,
0212                           80,
0213                           fBeamCenter - fBeamInterval,
0214                           fBeamCenter + fBeamInterval);
0215
0216   // Extrapolated hits
0217   fExtrapR = new TH1F("fExtrapR", "R^{extrap. hit}", 300, 400, 1000);
0218   fExtrapRho = new TH1F("fExtrapRho", "#rho^{extrap. hit}", 300, 400, 1000);
0219   fExtrapEta = new TH1F("fExtrapEta", "#eta^{extrap. hit}", 100, -0.1, 0.1);
0220   fExtrapPhi =
0221       new TH1F("fExtrapPhi", "#varphi^{extrap. hit}", 160, fBeamCenter - fBeamInterval, fBeamCenter + fBeamInterval);
0222
0223   fExtrapRVsPhi = new TH2F("fExtrapRVsPhi",
0224                            "R vs. #varphi for extrap. hits",
0225                            60,
0226                            400,
0227                            1000,
0228                            80,
0229                            fBeamCenter - fBeamInterval,
0230                            fBeamCenter + fBeamInterval);
0231
0232   // Distances
0233   fDeltaRo = new TH1F("fDeltaRo", "#Delta(#rho^{sim}, #rho^{extrap})", 100, 0, 200);
0234   fDeltaEta = new TH1F("fDeltaEta", "#Delta(#eta^{sim}, #eta^{extrap})", 100, -1, 1);
0235   fDeltaPhi = new TH1F("fDeltaPhi", "#Delta(#varphi^{sim}, #varphi^{extrap})", 120, -30, 30);
0236
0237   // Studies on Phi
0238   fStationPosPhi = new TH1F("fStationPosPhi", "Station with positive #varphi", 6, -0.5, 5.5);
0239   fSectorPosPhi = new TH1F("fSectorPosPhi", "Sector with positive #varphi", 6, -0.5, 5.5);
0240   fSLayerPosPhi = new TH1F("fSLayerPosPhi", "Superlayer with positive #varphi", 6, -0.5, 5.5);
0241   fLayerPosPhi = new TH1F("fLayerPosPhi", "Layer with positive #varphi", 6, -0.5, 5.5);
0242   fStationNegPhi = new TH1F("fStationNegPhi", "Station with negative #varphi", 6, -0.5, 5.5);
0243   fSectorNegPhi = new TH1F("fSectorNegPhi", "Sector with negative #varphi", 6, -0.5, 5.5);
0244   fSLayerNegPhi = new TH1F("fSLayerNegPhi", "Superlayer with negative #varphi", 6, -0.5, 5.5);
0245   fLayerNegPhi = new TH1F("fLayerNegPhi", "Layer with negative #varphi", 6, -0.5, 5.5);
0246
0247   //
0248   ///////////////////////////////////////////////////////////////////////////////////
0249 }
0250
0251 void Geant4ePropagatorAnalyzer::beginJob() { LogDebug("Geant4e") << "Nothing done in beginJob..."; }
0252
0253 void Geant4ePropagatorAnalyzer::endJob() {
0254   fDistanceSHLayer->Write();
0255   fDistance->Write();
0256   for (unsigned int i = 0; i < 4; i++)
0257     fDistanceSt[i]->Write();
0258
0259   fHitR->Write();
0260   fHitRho->Write();
0261   fHitEta->Write();
0262   fHitPhi->Write();
0263   fHitPhi1L->Write();
0264   fHitRVsPhi->Write();
0265
0266   fLayerR->Write();
0267   fLayerRho->Write();
0268   fLayerEta->Write();
0269   fLayerPhi->Write();
0270   fLayerRVsPhi->Write();
0271
0272   fExtrapR->Write();
0273   fExtrapRho->Write();
0274   fExtrapEta->Write();
0275   fExtrapPhi->Write();
0276   fExtrapRVsPhi->Write();
0277
0278   fDeltaRo->Write();
0279   fDeltaEta->Write();
0280   fDeltaPhi->Write();
0281
0282   // Phi studies
0283   fStationPosPhi->Write();
0284   fSectorPosPhi->Write();
0285   fSLayerPosPhi->Write();
0286   fLayerPosPhi->Write();
0287   fStationNegPhi->Write();
0288   fSectorNegPhi->Write();
0289   fSLayerNegPhi->Write();
0290   fLayerNegPhi->Write();
0291
0292   theRootFile->Close();
0293   //  TimingReport::current()->dump(std::cout);
0294 }
0295
0296 void Geant4ePropagatorAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0297   using namespace edm;
0298
0299   LogDebug("Geant4e") << "Starting analyze...";
0300
0301   ///////////////////////////////////////
0302   // Construct Magnetic Field
0303   const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
0304   if (bField.isValid())
0305     LogDebug("Geant4e") << "G4e -- Magnetic field is valid. Value in (0,0,0): "
0306                         << bField->inTesla(GlobalPoint(0, 0, 0)).mag() << " Tesla";
0307   else
0308     LogError("Geant4e") << "G4e -- NO valid Magnetic field";
0309
0310   ///////////////////////////////////////
0311   // Build geometry
0312
0313   //- DT...
0314   theDTGeomESH = iSetup.getHandle(dtGeomToken);
0315   LogDebug("Geant4e") << "Got DTGeometry";
0316
0317   //- CSC...
0318   theCSCGeomESH = iSetup.getHandle(cscGeomToken);
0319   LogDebug("Geant4e") << "Got CSCGeometry";
0320
0321   //- RPC...
0322   theRPCGeomESH = iSetup.getHandle(rpcGeomToken);
0323   LogDebug("Geant4e") << "Got RPCGeometry";
0324
0325   ///////////////////////////////////////
0326   // Run/Event information
0327   theRun = (int)iEvent.id().run();
0328   theEvent = (int)iEvent.id().event();
0329   LogDebug("Geant4e") << "G4e -- Begin for run/event ==" << theRun << "/" << theEvent
0330                       << " ---------------------------------";
0331
0332   ///////////////////////////////////////
0333   // Initialise the propagator
0334   if (!thePropagator)
0335     thePropagator = new Geant4ePropagator(&*bField);
0336
0337   if (thePropagator)
0338     LogDebug("Geant4e") << "Propagator built!";
0339   else
0340     LogError("Geant4e") << "Could not build propagator!";
0341
0342   ///////////////////////////////////////
0343   // Get the sim tracks & vertices
0344   Handle<SimTrackContainer> simTracks = iEvent.getHandle(simTrackToken_);
0345   if (!simTracks.isValid()) {
0346     LogWarning("Geant4e") << "No tracks found" << std::endl;
0347     return;
0348   }
0349   LogDebug("Geant4e") << "G4e -- Got simTracks of size " << simTracks->size();
0350
0351   Handle<SimVertexContainer> simVertices = iEvent.getHandle(simVertexToken_);
0352   if (!simVertices.isValid()) {
0353     LogWarning("Geant4e") << "No vertices found" << std::endl;
0354     return;
0355   }
0356   LogDebug("Geant4e") << "Got simVertices of size " << simVertices->size();
0357
0358   ///////////////////////////////////////
0359   // Get the sim hits for the different muon parts
0360   Handle<PSimHitContainer> simHitsDT = iEvent.getHandle(simHitsDTToken_);
0361   if (!simHitsDT.isValid()) {
0362     LogWarning("Geant4e") << "No hits found" << std::endl;
0363     return;
0364   }
0365   LogDebug("Geant4e") << "Got MuonDTHits of size " << simHitsDT->size();
0366
0367   Handle<PSimHitContainer> simHitsCSC = iEvent.getHandle(simHitsCSCToken_);
0368   if (!simHitsCSC.isValid()) {
0369     LogWarning("Geant4e") << "No hits found" << std::endl;
0370     return;
0371   }
0372   LogDebug("Geant4e") << "Got MuonCSCHits of size " << simHitsCSC->size();
0373
0374   Handle<PSimHitContainer> simHitsRPC = iEvent.getHandle(simHitsRPCToken_);
0375   if (!simHitsRPC.isValid()) {
0376     LogWarning("Geant4e") << "No hits found" << std::endl;
0377     return;
0378   }
0379   LogDebug("Geant4e") << "Got MuonRPCHits of size " << simHitsRPC->size();
0380
0381   ///////////////////////////////////////
0382   // Iterate over sim tracks to build the FreeTrajectoryState for
0383   // for the initial position.
0384   // DEBUG
0385   unsigned int counter = 0;
0386   // DEBUG
0387   for (SimTrackContainer::const_iterator simTracksIt = simTracks->begin(); simTracksIt != simTracks->end();
0388        simTracksIt++) {
0389     // DEBUG
0390     counter++;
0391     LogDebug("Geant4e") << "G4e -- Iterating over " << counter << "rd track. Number: " << simTracksIt->genpartIndex()
0392                         << "------------------";
0393     // DEBUG
0394
0395     int simTrackPDG = simTracksIt->type();
0396     if (abs(simTrackPDG) != 13) {
0397       continue;
0398     }
0399     LogDebug("Geant4e") << "G4e -- Track PDG " << simTrackPDG;
0400
0401     //- Timing
0402     //    TimeMe tProp("Geant4ePropagatorAnalyzer::analyze::propagate");
0403
0404     //- Check if the track corresponds to a muon
0405     int trkPDG = simTracksIt->type();
0406     if (abs(trkPDG) != 13) {
0407       LogDebug("Geant4e") << "Track is not a muon: " << trkPDG;
0408       continue;
0409     } else
0410       LogDebug("Geant4e") << "Found a muon track " << trkPDG;
0411
0412     //- Get momentum, but only use tracks with P > 2 GeV
0413     GlobalVector p3T = TrackPropagation::hep3VectorToGlobalVector(
0414         CLHEP::Hep3Vector(simTracksIt->momentum().x(), simTracksIt->momentum().y(), simTracksIt->momentum().z()));
0415     if (p3T.perp() < 2.) {
0416       LogDebug("Geant4e") << "Track PT is too low: " << p3T.perp();
0417       continue;
0418     } else {
0419       LogDebug("Geant4e") << "*** Phi (rad): " << p3T.phi() << " - Phi(deg)" << p3T.phi().degrees();
0420       LogDebug("Geant4e") << "Track PT is enough.";
0421       LogDebug("Geant4e") << "Track P.: " << p3T << "\nTrack P.: PT=" << p3T.perp() << "\tEta=" << p3T.eta()
0422                           << "\tPhi=" << p3T.phi().degrees() << "--> Rad: Phi=" << p3T.phi();
0423     }
0424
0425     //- Vertex fixes the starting point
0426     int vtxInd = simTracksIt->vertIndex();
0427     GlobalPoint r3T(0., 0., 0.);
0428     if (vtxInd < 0)
0429       LogDebug("Geant4e") << "Track with no vertex, defaulting to (0,0,0)";
0430     else
0431       // seems to be stored in mm --> convert to cm
0432       r3T = TrackPropagation::hep3VectorToGlobalPoint(CLHEP::Hep3Vector((*simVertices)[vtxInd].position().x(),
0433                                                                         (*simVertices)[vtxInd].position().y(),
0434                                                                         (*simVertices)[vtxInd].position().z()));
0435
0436     LogDebug("Geant4e") << "Init point: " << r3T << "\nInit point Ro=" << r3T.perp() << "\tEta=" << r3T.eta()
0437                         << "\tPhi=" << r3T.phi().degrees();
0438
0439     //- Charge
0440     int charge = trkPDG > 0 ? -1 : 1;
0441     LogDebug("Geant4e") << "Track charge = " << charge;
0442
0443     //- Initial covariance matrix is unity 10-6
0444     CurvilinearTrajectoryError covT;
0445     // covT *= 1E-6;
0446
0447     //- Build FreeTrajectoryState
0448     GlobalTrajectoryParameters trackPars(r3T, p3T, charge, &*bField);
0449     FreeTrajectoryState ftsTrack(trackPars, covT);
0450
0451     //- Get index of generated particle. Used further down
0452     unsigned int trkInd = simTracksIt->genpartIndex();
0453
0454     ////////////////////////////////////////////////
0455     //- Iterate over Sim Hits in DT and check propagation
0456     iterateOverHits(simHitsDT, DT, trkInd, ftsTrack);
0457     ////////////////////////////////////////////////
0458     //- Iterate over Sim Hits in RPC and check propagation
0459     // iterateOverHits(simHitsRPC, RPC, trkInd, ftsTrack);
0460     ////////////////////////////////////////////////
0461     //- Iterate over Sim Hits in CSC and check propagation
0462     // iterateOverHits(simHitsCSC, CSC, trkInd, ftsTrack);
0463
0464   }  // <-- for over sim tracks
0465 }
0466
0467 void Geant4ePropagatorAnalyzer::iterateOverHits(edm::Handle<edm::PSimHitContainer> simHits,
0468                                                 testMuChamberType muonChamberType,
0469                                                 unsigned int trkIndex,
0470                                                 const FreeTrajectoryState &ftsTrack) {
0471   using namespace edm;
0472
0473   if (muonChamberType == DT)
0474     LogDebug("Geant4e") << "G4e -- Iterating over DT hits...";
0475   else if (muonChamberType == RPC)
0476     LogDebug("Geant4e") << "G4e -- Iterating over RPC hits...";
0477   else if (muonChamberType == CSC)
0478     LogDebug("Geant4e") << "G4e -- Iterating over CSC hits...";
0479
0480   for (PSimHitContainer::const_iterator simHitIt = simHits->begin(); simHitIt != simHits->end(); simHitIt++) {
0481     ///////////////
0482     // Skip if this hit does not belong to the track
0483     if (simHitIt->trackId() != trkIndex) {
0484       LogDebug("Geant4e") << "Hit (in tr " << simHitIt->trackId() << ") does not belong to track " << trkIndex;
0485       continue;
0486     }
0487
0488     LogDebug("Geant4e") << "G4e -- Hit belongs to track " << trkIndex;
0489
0490     //////////////
0491     // Skip if it is not a muon (this is checked before also)
0492     int trkPDG = simHitIt->particleType();
0493     if (abs(trkPDG) != 13) {
0494       LogDebug("Geant4e") << "Associated track is not a muon: " << trkPDG;
0495       continue;
0496     }
0497     LogDebug("Geant4e") << "G4e -- Found a hit corresponding to a muon " << trkPDG;
0498
0499     //////////////////////////////////////////////////////////
0500     // Build the surface. This is different for DT, RPC, CSC
0501     // const GeomDetUnit* layer = 0;
0502     const GeomDet *layer = nullptr;
0503     // * DT
0504     DTWireId *wIdDT = nullptr;  // For DT holds the information about the chamber
0505     if (muonChamberType == DT) {
0506       wIdDT = new DTWireId(simHitIt->detUnitId());
0507       int station = wIdDT->station();
0508       LogDebug("Geant4e") << "G4e -- DT Chamber. Station: " << station << ". DetId: " << wIdDT->layerId();
0509
0510       if (fStudyStation != -1 && fStudyStation != station)
0511         continue;
0512
0513       layer = theDTGeomESH->layer(*wIdDT);
0514       if (layer == nullptr) {
0515         LogDebug("Geant4e") << "Failed to get detector unit";
0516         continue;
0517       }
0518     }
0519     // * RPC
0520     else if (muonChamberType == RPC) {
0521       RPCDetId rpcId(simHitIt->detUnitId());
0522       layer = theRPCGeomESH->idToDet(rpcId);
0523       if (layer == nullptr) {
0524         LogDebug("Geant4e") << "Failed to get detector unit";
0525         continue;
0526       }
0527     }
0528
0529     // * CSC
0530     else if (muonChamberType == CSC) {
0531       CSCDetId cscId(simHitIt->detUnitId());
0532       layer = theCSCGeomESH->idToDet(cscId);
0533       if (layer == nullptr) {
0534         LogDebug("Geant4e") << "Failed to get detector unit";
0535         continue;
0536       }
0537     }
0538
0539     const Surface &surf = layer->surface();
0540     if (layer->geographicalId() != DTLayerId(simHitIt->detUnitId()))
0541       LogError("Geant4e") << "ERROR: wrong DetId";
0542
0543     //==>DEBUG
0544     // const BoundPlane& bp = layer->surface();
0545     // const Bounds& bounds = bp.bounds();
0546     // LogDebug("Geant4e") << "Surface: length = " << bounds.length()
0547     //        << ", thickness = " << bounds.thickness()
0548     //      << ", width = " << bounds.width();
0549     //<==DEBUG
0550
0551     ////////////
0552     // Discard hits with very low momentum ???
0553     GlobalVector p3Hit = surf.toGlobal(simHitIt->momentumAtEntry());
0554     if (p3Hit.perp() < 0.5)
0555       continue;
0556     GlobalPoint posHit = surf.toGlobal(simHitIt->localPosition());
0557     Point3DBase<float, GlobalTag> surfpos = surf.position();
0558     LogDebug("Geant4e") << "G4e -- Layer position: " << surfpos << " cm"
0559                         << "\nG4e --                   Ro=" << surfpos.perp() << "\tEta=" << surfpos.eta()
0560                         << "\tPhi=" << surfpos.phi().degrees() << "deg";
0561     LogDebug("Geant4e") << "G4e -- Sim Hit position: " << posHit << " cm"
0562                         << "\nG4e --                   Ro=" << posHit.perp() << "\tEta=" << posHit.eta()
0563                         << "\tPhi=" << posHit.phi().degrees() << "deg"
0564                         << "\t localpos=" << simHitIt->localPosition();
0565
0566     const Plane *bp = dynamic_cast<const Plane *>(&surf);
0567     if (bp != nullptr) {
0568       Float_t distance = bp->localZ(posHit);
0569       LogDebug("Geant4e") << "\nG4e -- Distance from plane to sim hit: " << distance << "cm";
0570       fDistanceSHLayer->Fill(distance);
0571     } else {
0572       LogWarning("Geant4e") << "G4e -- Layer is not a Plane!!!";
0573       fDistanceSHLayer->Fill(-1);
0574     }
0575
0576     LogDebug("Geant4e") << "Sim Hit Momentum PT=" << p3Hit.perp() << "\tEta=" << p3Hit.eta()
0577                         << "\tPhi=" << p3Hit.phi().degrees() << "deg";
0578
0579     /////////////////////////////////////////
0580     // Propagate: Need to explicetely
0581     TrajectoryStateOnSurface tSOSDest = thePropagator->propagate(ftsTrack, surf);
0582
0583     /////////////////////
0584     // Get hit position and extrapolation position to compare
0585     GlobalPoint posExtrap = tSOSDest.freeState()->position();
0586     //     GlobalPoint posExtrap(posExtrap_prov.theta(),
0587     //            posExtrap_prov.phi()*12,
0588     //            posExtrap_prov.mag());
0589
0590     GlobalVector posDistance = posExtrap - posHit;
0591     float distance = posDistance.mag();
0592     float simhitphi = posHit.phi().degrees();
0593     float layerphi = surfpos.phi().degrees();
0594     float extrapphi = posExtrap.phi().degrees();
0595     LogDebug("Geant4e") << "G4e -- Difference between hit and final position: " << distance << " cm\n"
0596                         << "G4e -- Transversal difference between hit and final position: " << posDistance.perp()
0597                         << " cm.";
0598     LogDebug("Geant4e") << "G4e -- Extrapolated position:" << posExtrap << " cm\n"
0599                         << "G4e --       (Rho, eta, phi): (" << posExtrap.perp() << " cm, " << posExtrap.eta() << ", "
0600                         << extrapphi << " deg)";
0601     LogDebug("Geant4e") << "G4e --          Hit position: " << posHit << " cm\n"
0602                         << "G4e --       (Rho, eta, phi): (" << posHit.perp() << " cm, " << posHit.eta() << ", "
0603                         << simhitphi << " deg)";
0604
0605     fDistance->Fill(distance);
0606     fDistanceSt[wIdDT->station() - 1]->Fill(distance);
0607
0608     fHitR->Fill(posHit.mag());
0609     fHitRho->Fill(posHit.perp());
0610     fHitEta->Fill(posHit.eta());
0611     fHitPhi->Fill(simhitphi);
0612     if (posHit.perp() < 500)
0613       fHitPhi1L->Fill(simhitphi);
0614     fHitRVsPhi->Fill(posHit.mag(), simhitphi);
0615
0616     fLayerR->Fill(surfpos.mag());
0617     fLayerRho->Fill(surfpos.perp());
0618     fLayerEta->Fill(surfpos.eta());
0619     fLayerPhi->Fill(layerphi);
0620     fLayerRVsPhi->Fill(surfpos.mag(), layerphi);
0621
0622     fExtrapR->Fill(posExtrap.mag());
0623     fExtrapRho->Fill(posExtrap.perp());
0624     fExtrapEta->Fill(posExtrap.eta());
0625     fExtrapPhi->Fill(extrapphi);
0626     fExtrapRVsPhi->Fill(posExtrap.mag(), extrapphi);
0627
0628     fDeltaRo->Fill(posExtrap.perp() - posHit.perp());
0629     fDeltaEta->Fill(posExtrap.eta() - posHit.eta());
0630     fDeltaPhi->Fill(extrapphi - simhitphi);
0631
0632     if (wIdDT) {
0633       if (simhitphi > 0) {
0634         fStationPosPhi->Fill(wIdDT->station());
0635         fSectorPosPhi->Fill(wIdDT->sector());
0636         fSLayerPosPhi->Fill(wIdDT->superlayer());
0637         fLayerPosPhi->Fill(wIdDT->layer());
0638       } else {
0639         fStationNegPhi->Fill(wIdDT->station());
0640         fSectorNegPhi->Fill(wIdDT->sector());
0641         fSLayerNegPhi->Fill(wIdDT->superlayer());
0642         fLayerNegPhi->Fill(wIdDT->layer());
0643       }
0644       // Some cleaning...
0645       delete wIdDT;
0646     }
0647
0648   }  //<== For over simhits
0649 }
0650
0651 // define this as a plug-in
0652 DEFINE_FWK_MODULE(Geant4ePropagatorAnalyzer);