Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:22:11

0001 ////////////////////////////
0002 // Geometry Checklist     //
0003 // Maps with PSimHits     //
0004 //                        //
0005 // Nicola Pozzobon - 2012 //
0006 ////////////////////////////
0007 
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/PluginManager/interface/ModuleDef.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0015 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0016 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0017 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0019 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include <TH2D.h>
0022 
0023 //////////////////////////////
0024 //                          //
0025 //     CLASS DEFINITION     //
0026 //                          //
0027 //////////////////////////////
0028 
0029 class AnalyzerSimHitMaps : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0030   /// Public methods
0031 public:
0032   /// Constructor/destructor
0033   explicit AnalyzerSimHitMaps(const edm::ParameterSet& iConfig);
0034   ~AnalyzerSimHitMaps() override;
0035   // Typical methods used on Loops over events
0036   void beginJob() override;
0037   void endJob() override;
0038   void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0039   void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0040   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0041 
0042   /// Private methods and variables
0043 private:
0044   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> getTokenTrackerGeometry_;
0045 
0046   /// Global Position of SimHits
0047   TH2D* hSimHit_Barrel_XY;
0048   TH2D* hSimHit_Barrel_XY_Zoom;
0049   TH2D* hSimHit_Endcap_Fw_XY;
0050   TH2D* hSimHit_Endcap_Bw_XY;
0051   TH2D* hSimHit_RZ;
0052   TH2D* hSimHit_Endcap_Fw_RZ_Zoom;
0053   TH2D* hSimHit_Endcap_Bw_RZ_Zoom;
0054 
0055   std::map<std::string, TH2D*> m_hSimHit_Barrel_XY_Survey;
0056   std::map<std::string, TH2D*> m_hSimHit_RZ_Survey;
0057 };
0058 
0059 //////////////////////////////////
0060 //                              //
0061 //     CLASS IMPLEMENTATION     //
0062 //                              //
0063 //////////////////////////////////
0064 
0065 //////////////
0066 // CONSTRUCTOR
0067 AnalyzerSimHitMaps::AnalyzerSimHitMaps(edm::ParameterSet const& iConfig)
0068     : getTokenTrackerGeometry_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()) {
0069   /// Insert here what you need to initialize
0070   usesResource("TFileService");
0071 }
0072 
0073 /////////////
0074 // DESTRUCTOR
0075 AnalyzerSimHitMaps::~AnalyzerSimHitMaps() {
0076   /// Insert here what you need to delete
0077   /// when you close the class instance
0078 }
0079 
0080 //////////
0081 // END JOB
0082 void AnalyzerSimHitMaps::endJob() {
0083   /// Things to be done at the exit of the event Loop
0084   std::cerr << " AnalyzerSimHitMaps::endJob" << std::endl;
0085   /// End of things to be done at the exit from the event Loop
0086 }
0087 
0088 ////////////
0089 // BEGIN JOB
0090 void AnalyzerSimHitMaps::beginJob() {
0091   /// Initialize all slave variables
0092   /// mainly histogram ranges and resolution
0093   std::ostringstream histoName;
0094   std::ostringstream histoTitle;
0095 
0096   /// Things to be done before entering the event Loop
0097   std::cerr << " AnalyzerSimHitMaps::beginJob" << std::endl;
0098 
0099   /// Book histograms etc
0100   edm::Service<TFileService> fs;
0101 
0102   hSimHit_Barrel_XY = fs->make<TH2D>("hSimHit_Barrel_XY", "PSimHit Barrel y vs. x", 960, -120, 120, 960, -120, 120);
0103   hSimHit_Barrel_XY_Zoom =
0104       fs->make<TH2D>("hSimHit_Barrel_XY_Zoom", "PSimHit Barrel y vs. x", 960, 30, 60, 960, -15, 15);
0105   hSimHit_Endcap_Fw_XY =
0106       fs->make<TH2D>("hSimHit_Endcap_Fw_XY", "PSimHit Forward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0107   hSimHit_Endcap_Bw_XY =
0108       fs->make<TH2D>("hSimHit_Endcap_Bw_XY", "PSimHit Backward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0109   hSimHit_RZ = fs->make<TH2D>("hSimHit_RZ", "PSimHit #rho vs. z", 900, -300, 300, 480, 0, 120);
0110   hSimHit_Endcap_Fw_RZ_Zoom =
0111       fs->make<TH2D>("hSimHit_Endcap_Fw_RZ_Zoom", "PSimHit Forward Endcap #rho vs. z", 960, 140, 170, 960, 30, 60);
0112   hSimHit_Endcap_Bw_RZ_Zoom =
0113       fs->make<TH2D>("hSimHit_Endcap_Bw_RZ_Zoom", "PSimHit Backward Endcap #rho vs. z", 960, -170, -140, 960, 70, 100);
0114 
0115   hSimHit_Barrel_XY->Sumw2();
0116   hSimHit_Barrel_XY_Zoom->Sumw2();
0117   hSimHit_Endcap_Fw_XY->Sumw2();
0118   hSimHit_Endcap_Bw_XY->Sumw2();
0119   hSimHit_RZ->Sumw2();
0120   hSimHit_Endcap_Fw_RZ_Zoom->Sumw2();
0121   hSimHit_Endcap_Bw_RZ_Zoom->Sumw2();
0122 
0123   for (int ix = 0; ix < 11; ix++) {
0124     for (int iy = 10; iy > -1; iy--) {
0125       histoName.str("");
0126       histoTitle.str("");
0127       histoName << "hSimHit_Barrel_XY_Survey_" << -110 + ix * 20 << "x" << -110 + (1 + ix) * 20 << "_" << -110 + iy * 20
0128                 << "y" << -110 + (1 + iy) * 20;
0129       histoTitle << "PSimHit Barrel y (" << -110 + iy * 20 << ", " << -110 + (1 + iy) * 20 << ") cm vs. x ("
0130                  << -110 + ix * 20 << ", " << -110 + (1 + ix) * 20 << ") cm";
0131       TH2D* h = fs->make<TH2D>(histoName.str().c_str(),
0132                                histoTitle.str().c_str(),
0133                                800,
0134                                -110 + ix * 20,
0135                                -110 + (1 + ix) * 20,
0136                                800,
0137                                -110 + iy * 20,
0138                                -110 + (1 + iy) * 20);
0139       m_hSimHit_Barrel_XY_Survey.insert(std::pair<std::string, TH2D*>(histoName.str(), h));
0140     }
0141   }
0142 
0143   for (int iz = 0; iz < 27; iz++) {
0144     for (int ir = 5; ir >= 0; ir--) {
0145       histoName.str("");
0146       histoTitle.str("");
0147       histoName << "hSimHit_RZ_Survey_" << -10 + ir * 20 << "r" << -10 + (1 + ir) * 20 << "_" << -270 + iz * 20 << "z"
0148                 << -270 + (1 + iz) * 20;
0149       histoTitle << "PSimHit #rho (" << -10 + ir * 20 << ", " << -10 + (1 + ir) * 20 << ") cm vs. z (" << -270 + iz * 20
0150                  << ", " << -270 + (1 + iz) * 20 << ") cm";
0151       TH2D* h = fs->make<TH2D>(histoName.str().c_str(),
0152                                histoTitle.str().c_str(),
0153                                800,
0154                                -270 + iz * 20,
0155                                -270 + (1 + iz) * 20,
0156                                800,
0157                                -10 + ir * 20,
0158                                -10 + (1 + ir) * 20);
0159       m_hSimHit_RZ_Survey.insert(std::pair<std::string, TH2D*>(histoName.str(), h));
0160     }
0161   }
0162 
0163   /// End of things to be done before entering the event Loop
0164 }
0165 
0166 //////////
0167 // ANALYZE
0168 void AnalyzerSimHitMaps::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0169   /// Geometry handles etc
0170   edm::ESHandle<TrackerGeometry> geometryHandle = iSetup.getHandle(getTokenTrackerGeometry_);
0171   const TrackerGeometry* theGeometry = &(*geometryHandle);
0172 
0173   //////////////////
0174   // GET SIM HITS //
0175   //////////////////
0176   edm::Handle<edm::PSimHitContainer> simHitHandlePXBH;
0177   edm::Handle<edm::PSimHitContainer> simHitHandlePXBL;
0178   edm::Handle<edm::PSimHitContainer> simHitHandlePXFH;
0179   edm::Handle<edm::PSimHitContainer> simHitHandlePXFL;
0180   iEvent.getByLabel("g4SimHits", "TrackerHitsPixelBarrelHighTof", simHitHandlePXBH);
0181   iEvent.getByLabel("g4SimHits", "TrackerHitsPixelBarrelLowTof", simHitHandlePXBL);
0182   iEvent.getByLabel("g4SimHits", "TrackerHitsPixelEndcapHighTof", simHitHandlePXFH);
0183   iEvent.getByLabel("g4SimHits", "TrackerHitsPixelEndcapLowTof", simHitHandlePXFL);
0184 
0185   /// Loop over Sim Hits collections
0186   edm::PSimHitContainer::const_iterator iterSimHit;
0187   for (iterSimHit = simHitHandlePXBH->begin(); iterSimHit != simHitHandlePXBH->end(); iterSimHit++) {
0188     const GeomDetUnit* gDetUnit = theGeometry->idToDetUnit(iterSimHit->detUnitId());
0189     GlobalPoint shPos = gDetUnit->surface().toGlobal(iterSimHit->localPosition());
0190     std::ostringstream histoNameXY;
0191     std::ostringstream histoNameRZ;
0192     histoNameXY << "hSimHit_Barrel_XY_Survey_" << (20 * (floor((shPos.x() + 10) / 20)) - 10) << "x"
0193                 << (20 * (floor((shPos.x() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.y() + 10) / 20)) - 10)
0194                 << "y" << (20 * (floor((shPos.y() + 10) / 20)) + 10);
0195     histoNameRZ << "hSimHit_RZ_Survey_" << (20 * (floor((shPos.perp() + 10) / 20)) - 10) << "r"
0196                 << (20 * (floor((shPos.perp() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.z() + 10) / 20)) - 10)
0197                 << "z" << (20 * (floor((shPos.z() + 10) / 20)) + 10);
0198     hSimHit_RZ->Fill(shPos.z(), shPos.perp());
0199     if (m_hSimHit_RZ_Survey.find(histoNameRZ.str()) != m_hSimHit_RZ_Survey.end())
0200       m_hSimHit_RZ_Survey.find(histoNameRZ.str())->second->Fill(shPos.z(), shPos.perp());
0201     if (gDetUnit->type().isBarrel()) {
0202       hSimHit_Barrel_XY->Fill(shPos.x(), shPos.y());
0203       hSimHit_Barrel_XY_Zoom->Fill(shPos.x(), shPos.y());
0204       if (m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str()) != m_hSimHit_Barrel_XY_Survey.end())
0205         m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str())->second->Fill(shPos.x(), shPos.y());
0206     } else if (gDetUnit->type().isEndcap()) {
0207       if (shPos.z() > 0) {
0208         hSimHit_Endcap_Fw_XY->Fill(shPos.x(), shPos.y());
0209         hSimHit_Endcap_Fw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0210       } else {
0211         hSimHit_Endcap_Bw_XY->Fill(shPos.x(), shPos.y());
0212         hSimHit_Endcap_Bw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0213       }
0214     }
0215   }
0216 
0217   for (iterSimHit = simHitHandlePXBL->begin(); iterSimHit != simHitHandlePXBL->end(); iterSimHit++) {
0218     const GeomDetUnit* gDetUnit = theGeometry->idToDetUnit(iterSimHit->detUnitId());
0219     GlobalPoint shPos = gDetUnit->surface().toGlobal(iterSimHit->localPosition());
0220     std::ostringstream histoNameXY;
0221     std::ostringstream histoNameRZ;
0222     histoNameXY << "hSimHit_Barrel_XY_Survey_" << (20 * (floor((shPos.x() + 10) / 20)) - 10) << "x"
0223                 << (20 * (floor((shPos.x() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.y() + 10) / 20)) - 10)
0224                 << "y" << (20 * (floor((shPos.y() + 10) / 20)) + 10);
0225     histoNameRZ << "hSimHit_RZ_Survey_" << (20 * (floor((shPos.perp() + 10) / 20)) - 10) << "r"
0226                 << (20 * (floor((shPos.perp() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.z() + 10) / 20)) - 10)
0227                 << "z" << (20 * (floor((shPos.z() + 10) / 20)) + 10);
0228     hSimHit_RZ->Fill(shPos.z(), shPos.perp());
0229     if (m_hSimHit_RZ_Survey.find(histoNameRZ.str()) != m_hSimHit_RZ_Survey.end())
0230       m_hSimHit_RZ_Survey.find(histoNameRZ.str())->second->Fill(shPos.z(), shPos.perp());
0231     if (gDetUnit->type().isBarrel()) {
0232       hSimHit_Barrel_XY->Fill(shPos.x(), shPos.y());
0233       hSimHit_Barrel_XY_Zoom->Fill(shPos.x(), shPos.y());
0234       if (m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str()) != m_hSimHit_Barrel_XY_Survey.end())
0235         m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str())->second->Fill(shPos.x(), shPos.y());
0236     } else if (gDetUnit->type().isEndcap()) {
0237       if (shPos.z() > 0) {
0238         hSimHit_Endcap_Fw_XY->Fill(shPos.x(), shPos.y());
0239         hSimHit_Endcap_Fw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0240       } else {
0241         hSimHit_Endcap_Bw_XY->Fill(shPos.x(), shPos.y());
0242         hSimHit_Endcap_Bw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0243       }
0244     }
0245   }
0246 
0247   for (iterSimHit = simHitHandlePXFH->begin(); iterSimHit != simHitHandlePXFH->end(); iterSimHit++) {
0248     const GeomDetUnit* gDetUnit = theGeometry->idToDetUnit(iterSimHit->detUnitId());
0249     GlobalPoint shPos = gDetUnit->surface().toGlobal(iterSimHit->localPosition());
0250     std::ostringstream histoNameXY;
0251     std::ostringstream histoNameRZ;
0252     histoNameXY << "hSimHit_Barrel_XY_Survey_" << (20 * (floor((shPos.x() + 10) / 20)) - 10) << "x"
0253                 << (20 * (floor((shPos.x() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.y() + 10) / 20)) - 10)
0254                 << "y" << (20 * (floor((shPos.y() + 10) / 20)) + 10);
0255     histoNameRZ << "hSimHit_RZ_Survey_" << (20 * (floor((shPos.perp() + 10) / 20)) - 10) << "r"
0256                 << (20 * (floor((shPos.perp() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.z() + 10) / 20)) - 10)
0257                 << "z" << (20 * (floor((shPos.z() + 10) / 20)) + 10);
0258     hSimHit_RZ->Fill(shPos.z(), shPos.perp());
0259     if (m_hSimHit_RZ_Survey.find(histoNameRZ.str()) != m_hSimHit_RZ_Survey.end())
0260       m_hSimHit_RZ_Survey.find(histoNameRZ.str())->second->Fill(shPos.z(), shPos.perp());
0261     if (gDetUnit->type().isBarrel()) {
0262       hSimHit_Barrel_XY->Fill(shPos.x(), shPos.y());
0263       hSimHit_Barrel_XY_Zoom->Fill(shPos.x(), shPos.y());
0264       if (m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str()) != m_hSimHit_Barrel_XY_Survey.end())
0265         m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str())->second->Fill(shPos.x(), shPos.y());
0266     } else if (gDetUnit->type().isEndcap()) {
0267       if (shPos.z() > 0) {
0268         hSimHit_Endcap_Fw_XY->Fill(shPos.x(), shPos.y());
0269         hSimHit_Endcap_Fw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0270       } else {
0271         hSimHit_Endcap_Bw_XY->Fill(shPos.x(), shPos.y());
0272         hSimHit_Endcap_Bw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0273       }
0274     }
0275   }
0276 
0277   for (iterSimHit = simHitHandlePXFL->begin(); iterSimHit != simHitHandlePXFL->end(); iterSimHit++) {
0278     const GeomDetUnit* gDetUnit = theGeometry->idToDetUnit(iterSimHit->detUnitId());
0279     GlobalPoint shPos = gDetUnit->surface().toGlobal(iterSimHit->localPosition());
0280     std::ostringstream histoNameXY;
0281     std::ostringstream histoNameRZ;
0282     histoNameXY << "hSimHit_Barrel_XY_Survey_" << (20 * (floor((shPos.x() + 10) / 20)) - 10) << "x"
0283                 << (20 * (floor((shPos.x() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.y() + 10) / 20)) - 10)
0284                 << "y" << (20 * (floor((shPos.y() + 10) / 20)) + 10);
0285     histoNameRZ << "hSimHit_RZ_Survey_" << (20 * (floor((shPos.perp() + 10) / 20)) - 10) << "r"
0286                 << (20 * (floor((shPos.perp() + 10) / 20)) + 10) << "_" << (20 * (floor((shPos.z() + 10) / 20)) - 10)
0287                 << "z" << (20 * (floor((shPos.z() + 10) / 20)) + 10);
0288     hSimHit_RZ->Fill(shPos.z(), shPos.perp());
0289     if (m_hSimHit_RZ_Survey.find(histoNameRZ.str()) != m_hSimHit_RZ_Survey.end())
0290       m_hSimHit_RZ_Survey.find(histoNameRZ.str())->second->Fill(shPos.z(), shPos.perp());
0291     if (gDetUnit->type().isBarrel()) {
0292       hSimHit_Barrel_XY->Fill(shPos.x(), shPos.y());
0293       hSimHit_Barrel_XY_Zoom->Fill(shPos.x(), shPos.y());
0294       if (m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str()) != m_hSimHit_Barrel_XY_Survey.end())
0295         m_hSimHit_Barrel_XY_Survey.find(histoNameXY.str())->second->Fill(shPos.x(), shPos.y());
0296     } else if (gDetUnit->type().isEndcap()) {
0297       if (shPos.z() > 0) {
0298         hSimHit_Endcap_Fw_XY->Fill(shPos.x(), shPos.y());
0299         hSimHit_Endcap_Fw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0300       } else {
0301         hSimHit_Endcap_Bw_XY->Fill(shPos.x(), shPos.y());
0302         hSimHit_Endcap_Bw_RZ_Zoom->Fill(shPos.z(), shPos.perp());
0303       }
0304     }
0305   }
0306 
0307 }  /// End of analyze()
0308 
0309 ///////////////////////////
0310 // DEFINE THIS AS A PLUG-IN
0311 DEFINE_FWK_MODULE(AnalyzerSimHitMaps);