Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //- 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 {}
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   // 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   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   // debug_ = iConfig.getParameter<bool>("debug");
0151   fStudyStation = iConfig.getParameter<int>("StudyStation");
0152 
0153   fBeamCenter = iConfig.getParameter<double>("BeamCenter");
0154   fBeamInterval = iConfig.getParameter<double>("BeamInterval");
0155 
0156   ///////////////////////////////////////////////////////////////////////////////////
0157   // Histograms
0158   //
0159 
0160   // Output ROOT file
0161   theRootFile = new TFile(iConfig.getParameter<std::string>("RootFile").c_str(), "recreate");
0162 
0163   // Distance between Sim Hit and associated Layer
0164   fDistanceSHLayer = new TH1F("fDistanceSHLayer", "Distance(sim hit - layer)", 200, -1, 1);
0165 
0166   // Distance between simhit and extrapolation
0167   fDistance = new TH1F("fDistance", "Distance(sim hit - extrap)", 150, 0, 300);
0168   // Distance between simhit and extrapolation
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   // Simulated hits
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   // Position of layers
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   // Extrapolated hits
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   // Distances
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   // Studies on Phi
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   // Phi studies
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   //  TimingReport::current()->dump(std::cout);
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   // Construct Magnetic Field
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   // Build geometry
0301 
0302   //- DT...
0303   theDTGeomESH = iSetup.getHandle(dtGeomToken);
0304   LogDebug("Geant4e") << "Got DTGeometry";
0305 
0306   //- CSC...
0307   theCSCGeomESH = iSetup.getHandle(cscGeomToken);
0308   LogDebug("Geant4e") << "Got CSCGeometry";
0309 
0310   //- RPC...
0311   theRPCGeomESH = iSetup.getHandle(rpcGeomToken);
0312   LogDebug("Geant4e") << "Got RPCGeometry";
0313 
0314   ///////////////////////////////////////
0315   // Run/Event information
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   // Initialise the propagator
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   // Get the sim tracks & vertices
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   // Get the sim hits for the different muon parts
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   // Iterate over sim tracks to build the FreeTrajectoryState for
0377   // for the initial position.
0378   // DEBUG
0379   unsigned int counter = 0;
0380   // DEBUG
0381   for (SimTrackContainer::const_iterator simTracksIt = simTracks->begin(); simTracksIt != simTracks->end();
0382        simTracksIt++) {
0383     // DEBUG
0384     counter++;
0385     LogDebug("Geant4e") << "G4e -- Iterating over " << counter << "rd track. Number: " << simTracksIt->genpartIndex()
0386                         << "------------------";
0387     // DEBUG
0388 
0389     int simTrackPDG = simTracksIt->type();
0390     if (abs(simTrackPDG) != 13) {
0391       continue;
0392     }
0393     LogDebug("Geant4e") << "G4e -- Track PDG " << simTrackPDG;
0394 
0395     //- Timing
0396     //    TimeMe tProp("Geant4ePropagatorAnalyzer::analyze::propagate");
0397 
0398     //- Check if the track corresponds to a muon
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     //- Get momentum, but only use tracks with P > 2 GeV
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     //- Vertex fixes the starting point
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       // seems to be stored in mm --> convert to cm
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     //- Charge
0434     int charge = trkPDG > 0 ? -1 : 1;
0435     LogDebug("Geant4e") << "Track charge = " << charge;
0436 
0437     //- Initial covariance matrix is unity 10-6
0438     CurvilinearTrajectoryError covT;
0439     // covT *= 1E-6;
0440 
0441     //- Build FreeTrajectoryState
0442     GlobalTrajectoryParameters trackPars(r3T, p3T, charge, &*bField);
0443     FreeTrajectoryState ftsTrack(trackPars, covT);
0444 
0445     //- Get index of generated particle. Used further down
0446     unsigned int trkInd = simTracksIt->genpartIndex();
0447 
0448     ////////////////////////////////////////////////
0449     //- Iterate over Sim Hits in DT and check propagation
0450     iterateOverHits(simHitsDT, DT, trkInd, ftsTrack);
0451     ////////////////////////////////////////////////
0452     //- Iterate over Sim Hits in RPC and check propagation
0453     // iterateOverHits(simHitsRPC, RPC, trkInd, ftsTrack);
0454     ////////////////////////////////////////////////
0455     //- Iterate over Sim Hits in CSC and check propagation
0456     // iterateOverHits(simHitsCSC, CSC, trkInd, ftsTrack);
0457 
0458   }  // <-- for over sim tracks
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     // Skip if this hit does not belong to the track
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     // Skip if it is not a muon (this is checked before also)
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     // Build the surface. This is different for DT, RPC, CSC
0495     // const GeomDetUnit* layer = 0;
0496     const GeomDet *layer = nullptr;
0497     // * DT
0498     DTWireId *wIdDT = nullptr;  // For DT holds the information about the chamber
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     // * RPC
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     // * CSC
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     //==>DEBUG
0538     // const BoundPlane& bp = layer->surface();
0539     // const Bounds& bounds = bp.bounds();
0540     // LogDebug("Geant4e") << "Surface: length = " << bounds.length()
0541     //        << ", thickness = " << bounds.thickness()
0542     //      << ", width = " << bounds.width();
0543     //<==DEBUG
0544 
0545     ////////////
0546     // Discard hits with very low momentum ???
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     // Propagate: Need to explicetely
0575     TrajectoryStateOnSurface tSOSDest = thePropagator->propagate(ftsTrack, surf);
0576 
0577     /////////////////////
0578     // Get hit position and extrapolation position to compare
0579     GlobalPoint posExtrap = tSOSDest.freeState()->position();
0580     //     GlobalPoint posExtrap(posExtrap_prov.theta(),
0581     //            posExtrap_prov.phi()*12,
0582     //            posExtrap_prov.mag());
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       // Some cleaning...
0639       delete wIdDT;
0640     }
0641 
0642   }  //<== For over simhits
0643 }
0644 
0645 // define this as a plug-in
0646 DEFINE_FWK_MODULE(Geant4ePropagatorAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!