Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:02

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