Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-27 03:18:10

0001 #include "Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h"
0002 
0003 #include "TFile.h"
0004 #include "TTree.h"
0005 #include "TBranch.h"
0006 #include "TH1F.h"
0007 
0008 #include <iostream>
0009 #include <string>
0010 
0011 using namespace edm;
0012 using namespace std;
0013 
0014 MuonSimHitsValidAnalyzer::MuonSimHitsValidAnalyzer(const edm::ParameterSet& iPSet)
0015     : fName(""),
0016       verbosity(0),
0017       label(""),
0018       getAllProvenances(false),
0019       printProvenanceInfo(false),
0020       nRawGenPart(0),
0021       count(0)
0022 
0023 {
0024   /// get information from parameter set
0025   fName = iPSet.getUntrackedParameter<std::string>("Name");
0026   verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0027   label = iPSet.getParameter<std::string>("Label");
0028   edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0029   getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0030   printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0031 
0032   nRawGenPart = 0;
0033 
0034   /// get labels for input tags
0035   DTHitsToken_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("DTHitsSrc"));
0036 
0037   /// print out Parameter Set information being used
0038 
0039   if (verbosity) {
0040     Labels l;
0041     labelsForToken(DTHitsToken_, l);
0042     edm::LogInfo("MuonSimHitsValidAnalyzer::MuonSimHitsValidAnalyzer")
0043         << "\n===============================\n"
0044         << "Initialized as EDAnalyzer with parameter values:\n"
0045         << "    Name      = " << fName << "\n"
0046         << "    Verbosity = " << verbosity << "\n"
0047         << "    Label     = " << label << "\n"
0048         << "    GetProv   = " << getAllProvenances << "\n"
0049         << "    PrintProv = " << printProvenanceInfo
0050         << "\n"
0051         //    << "    CSCHitsSrc=  " <<CSCHitsSrc_.label()
0052         //    << ":" << CSCHitsSrc_.instance() << "\n"
0053         << "    DTHitsSrc =  " << l.module << ":" << l.productInstance
0054         << "\n"
0055         //     << "    RPCHitsSrc=  " <<RPCHitsSrc_.label()
0056         //     << ":" << RPCHitsSrc_.instance() << "\n"
0057         << "===============================\n";
0058   }
0059 
0060   pow6 = 1000000.0;
0061   mom4 = 0.;
0062   mom1 = 0;
0063   costeta = 0.;
0064   radius = 0;
0065   sinteta = 0.;
0066   globposx = 0.;
0067   globposy = 0;
0068   nummu_DT = 0;
0069   nummu_CSC = 0;
0070   nummu_RPC = 0;
0071   geomToken_ = esConsumes<DTGeometry, MuonGeometryRecord>();
0072 }
0073 
0074 MuonSimHitsValidAnalyzer::~MuonSimHitsValidAnalyzer() {}
0075 
0076 void MuonSimHitsValidAnalyzer::bookHistograms(DQMStore::IBooker& iBooker,
0077                                               edm::Run const& iRun,
0078                                               edm::EventSetup const& /* iSetup */) {
0079   meAllDTHits = nullptr;
0080   meMuDTHits = nullptr;
0081   meToF = nullptr;
0082   meEnergyLoss = nullptr;
0083   meMomentumMB1 = nullptr;
0084   meMomentumMB4 = nullptr;
0085   meLossMomIron = nullptr;
0086   meLocalXvsZ = nullptr;
0087   meLocalXvsY = nullptr;
0088   meGlobalXvsZ = nullptr;
0089   meGlobalXvsY = nullptr;
0090   meGlobalXvsZWm2 = nullptr;
0091   meGlobalXvsZWm1 = nullptr;
0092   meGlobalXvsZW0 = nullptr;
0093   meGlobalXvsZWp1 = nullptr;
0094   meGlobalXvsZWp2 = nullptr;
0095   meGlobalXvsYWm2 = nullptr;
0096   meGlobalXvsYWm1 = nullptr;
0097   meGlobalXvsYW0 = nullptr;
0098   meGlobalXvsYWp1 = nullptr;
0099   meGlobalXvsYWp2 = nullptr;
0100   meWheelOccup = nullptr;
0101   meStationOccup = nullptr;
0102   meSectorOccup = nullptr;
0103   meSuperLOccup = nullptr;
0104   meLayerOccup = nullptr;
0105   meWireOccup = nullptr;
0106   mePathMuon = nullptr;
0107   meChamberOccup = nullptr;
0108   meHitRadius = nullptr;
0109   meCosTheta = nullptr;
0110   meGlobalEta = nullptr;
0111   meGlobalPhi = nullptr;
0112 
0113   Char_t histo_n[100];
0114   Char_t histo_t[100];
0115 
0116   iBooker.setCurrentFolder("MuonDTHitsV/DTHitsValidationTask");
0117 
0118   sprintf(histo_n, "Number_of_all_DT_hits");
0119   sprintf(histo_t, "Number_of_all_DT_hits");
0120   meAllDTHits = iBooker.book1D(histo_n, histo_t, 200, 1.0, 201.0);
0121 
0122   sprintf(histo_n, "Number_of_muon_DT_hits");
0123   sprintf(histo_t, "Number_of_muon_DT_hits");
0124   meMuDTHits = iBooker.book1D(histo_n, histo_t, 150, 1.0, 151.0);
0125 
0126   sprintf(histo_n, "Tof_of_hits");
0127   sprintf(histo_t, "Tof_of_hits");
0128   meToF = iBooker.book1D(histo_n, histo_t, 100, -0.5, 50.);
0129 
0130   sprintf(histo_n, "DT_energy_loss_keV");
0131   sprintf(histo_t, "DT_energy_loss_keV");
0132   meEnergyLoss = iBooker.book1D(histo_n, histo_t, 100, 0.0, 10.0);
0133 
0134   sprintf(histo_n, "Momentum_at_MB1");
0135   sprintf(histo_t, "Momentum_at_MB1");
0136   meMomentumMB1 = iBooker.book1D(histo_n, histo_t, 100, 10.0, 200.0);
0137 
0138   sprintf(histo_n, "Momentum_at_MB4");
0139   sprintf(histo_t, "Momentum_at_MB4");
0140   meMomentumMB4 = iBooker.book1D(histo_n, histo_t, 100, 10.0, 200.0);
0141 
0142   sprintf(histo_n, "Loss_of_muon_Momentum_in_Iron");
0143   sprintf(histo_t, "Loss_of_muon_Momentum_in_Iron");
0144   meLossMomIron = iBooker.book1D(histo_n, histo_t, 80, 0.0, 40.0);
0145 
0146   sprintf(histo_n, "Local_x-coord_vs_local_z-coord_of_muon_hit");
0147   sprintf(histo_t, "Local_x-coord_vs_local_z-coord_of_muon_hit");
0148   meLocalXvsZ = iBooker.book2D(histo_n, histo_t, 100, -150., 150., 100, -0.8, 0.8);
0149 
0150   sprintf(histo_n, "local_x-coord_vs_local_y-coord_of_muon_hit");
0151   sprintf(histo_t, "local_x-coord_vs_local_y-coord_of_muon_hit");
0152   meLocalXvsY = iBooker.book2D(histo_n, histo_t, 100, -150., 150., 100, -150., 150.);
0153 
0154   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit");
0155   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit");
0156   meGlobalXvsZ = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0157 
0158   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit");
0159   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit");
0160   meGlobalXvsY = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0161 
0162   // New histos
0163 
0164   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-2");
0165   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-2");
0166   meGlobalXvsZWm2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0167 
0168   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-2");
0169   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-2");
0170   meGlobalXvsYWm2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0171 
0172   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-1");
0173   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-1");
0174   meGlobalXvsZWm1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0175 
0176   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-1");
0177   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-1");
0178   meGlobalXvsYWm1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0179 
0180   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w0");
0181   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w0");
0182   meGlobalXvsZW0 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0183 
0184   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w0");
0185   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w0");
0186   meGlobalXvsYW0 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0187 
0188   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w1");
0189   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w1");
0190   meGlobalXvsZWp1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0191 
0192   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w1");
0193   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w1");
0194   meGlobalXvsYWp1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0195 
0196   sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w2");
0197   sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w2");
0198   meGlobalXvsZWp2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0199 
0200   sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w2");
0201   sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w2");
0202   meGlobalXvsYWp2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);
0203 
0204   //
0205 
0206   sprintf(histo_n, "Wheel_occupancy");
0207   sprintf(histo_t, "Wheel_occupancy");
0208   meWheelOccup = iBooker.book1D(histo_n, histo_t, 10, -5.0, 5.0);
0209 
0210   sprintf(histo_n, "Station_occupancy");
0211   sprintf(histo_t, "Station_occupancy");
0212   meStationOccup = iBooker.book1D(histo_n, histo_t, 6, 0., 6.0);
0213 
0214   sprintf(histo_n, "Sector_occupancy");
0215   sprintf(histo_t, "Sector_occupancy");
0216   meSectorOccup = iBooker.book1D(histo_n, histo_t, 20, 0., 20.);
0217 
0218   sprintf(histo_n, "SuperLayer_occupancy");
0219   sprintf(histo_t, "SuperLayer_occupancy");
0220   meSuperLOccup = iBooker.book1D(histo_n, histo_t, 5, 0., 5.);
0221 
0222   sprintf(histo_n, "Layer_occupancy");
0223   sprintf(histo_t, "Layer_occupancy");
0224   meLayerOccup = iBooker.book1D(histo_n, histo_t, 6, 0., 6.);
0225 
0226   sprintf(histo_n, "Wire_occupancy");
0227   sprintf(histo_t, "Wire_occupancy");
0228   meWireOccup = iBooker.book1D(histo_n, histo_t, 100, 0., 100.);
0229 
0230   sprintf(histo_n, "path_followed_by_muon");
0231   sprintf(histo_t, "path_followed_by_muon");
0232   mePathMuon = iBooker.book1D(histo_n, histo_t, 160, 0., 160.);
0233 
0234   sprintf(histo_n, "chamber_occupancy");
0235   sprintf(histo_t, "chamber_occupancy");
0236   meChamberOccup = iBooker.book1D(histo_n, histo_t, 251, 0., 251.);
0237 
0238   sprintf(histo_n, "radius_of_hit");
0239   sprintf(histo_t, "radius_of_hit");
0240   meHitRadius = iBooker.book1D(histo_n, histo_t, 100, 0., 1200.);
0241 
0242   sprintf(histo_n, "costheta_of_hit");
0243   sprintf(histo_t, "costheta_of_hit");
0244   meCosTheta = iBooker.book1D(histo_n, histo_t, 100, -1., 1.);
0245 
0246   sprintf(histo_n, "global_eta_of_hit");
0247   sprintf(histo_t, "global_eta_of_hit");
0248   meGlobalEta = iBooker.book1D(histo_n, histo_t, 60, -2.7, 2.7);
0249 
0250   sprintf(histo_n, "global_phi_of_hit");
0251   sprintf(histo_t, "global_phi_of_hit");
0252   meGlobalPhi = iBooker.book1D(histo_n, histo_t, 60, -3.14, 3.14);
0253 }
0254 
0255 void MuonSimHitsValidAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0256   /// keep track of number of events processed
0257   ++count;
0258 
0259   /// get event id information
0260   edm::RunNumber_t nrun = iEvent.id().run();
0261   edm::EventNumber_t nevt = iEvent.id().event();
0262 
0263   if (verbosity > 0) {
0264     edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Processing run " << nrun << ", event " << nevt;
0265   }
0266 
0267   /// look at information available in the event
0268   if (getAllProvenances) {
0269     std::vector<const edm::StableProvenance*> AllProv;
0270     iEvent.getAllStableProvenance(AllProv);
0271 
0272     if (verbosity > 0)
0273       edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Number of Provenances = " << AllProv.size();
0274 
0275     if (printProvenanceInfo && (verbosity > 0)) {
0276       TString eventout("\nProvenance info:\n");
0277 
0278       for (unsigned int i = 0; i < AllProv.size(); ++i) {
0279         eventout += "\n       ******************************";
0280         eventout += "\n       Module       : ";
0281         eventout += AllProv[i]->moduleLabel();
0282         eventout += "\n       ProductID    : ";
0283         eventout += AllProv[i]->productID().id();
0284         eventout += "\n       ClassName    : ";
0285         eventout += AllProv[i]->className();
0286         eventout += "\n       InstanceName : ";
0287         eventout += AllProv[i]->productInstanceName();
0288         eventout += "\n       BranchName   : ";
0289         eventout += AllProv[i]->branchName();
0290       }
0291       eventout += "       ******************************\n";
0292       edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << eventout << "\n";
0293     }
0294   }
0295 
0296   fillDT(iEvent, iSetup);
0297 
0298   if (verbosity > 0)
0299     edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Done gathering data from event.";
0300 
0301   return;
0302 }
0303 
0304 void MuonSimHitsValidAnalyzer::fillDT(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0305   TString eventout;
0306   if (verbosity > 0)
0307     eventout = "\nGathering DT info:";
0308 
0309   /// iterator to access containers
0310   edm::PSimHitContainer::const_iterator itHit;
0311 
0312   /// access the DT
0313   edm::ESHandle<DTGeometry> hGeom = iSetup.getHandle(geomToken_);
0314   if (!hGeom.isValid()) {
0315     edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT")
0316         << "Unable to find MuonGeometryRecord for the DTGeometry in event!";
0317     return;
0318   }
0319   const DTGeometry& theDTMuon(*hGeom);
0320 
0321   /// get DT information
0322   edm::Handle<edm::PSimHitContainer> MuonDTContainer;
0323   iEvent.getByToken(DTHitsToken_, MuonDTContainer);
0324   if (!MuonDTContainer.isValid()) {
0325     edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT") << "Unable to find MuonDTHits in event!";
0326     return;
0327   }
0328 
0329   touch1 = 0;
0330   touch4 = 0;
0331   nummu_DT = 0;
0332 
0333   meAllDTHits->Fill(MuonDTContainer->size());
0334 
0335   /// cycle through container
0336   int i = 0, j = 0;
0337   for (itHit = MuonDTContainer->begin(); itHit != MuonDTContainer->end(); ++itHit) {
0338     ++i;
0339 
0340     /// create a DetId from the detUnitId
0341     DetId theDetUnitId(itHit->detUnitId());
0342     int detector = theDetUnitId.det();
0343     int subdetector = theDetUnitId.subdetId();
0344 
0345     /// check that expected detector is returned
0346     if ((detector == dMuon) && (subdetector == sdMuonDT)) {
0347       /// get the GeomDetUnit from the geometry using theDetUnitID
0348       const GeomDetUnit* theDet = theDTMuon.idToDetUnit(theDetUnitId);
0349 
0350       if (!theDet) {
0351         edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT") << "Unable to get GeomDetUnit from theDTMuon for hit " << i;
0352         continue;
0353       }
0354 
0355       ++j;
0356 
0357       /// get the Surface of the hit (knows how to go from local <-> global)
0358       const BoundPlane& bsurf = theDet->surface();
0359 
0360       /// gather necessary information
0361 
0362       if (abs(itHit->particleType()) == 13) {
0363         nummu_DT++;
0364         meToF->Fill(itHit->tof());
0365         meEnergyLoss->Fill(itHit->energyLoss() * pow6);
0366 
0367         iden = itHit->detUnitId();
0368 
0369         wheel = ((iden >> 15) & 0x7) - 3;
0370         station = ((iden >> 22) & 0x7);
0371         sector = ((iden >> 18) & 0xf);
0372         superlayer = ((iden >> 13) & 0x3);
0373         layer = ((iden >> 10) & 0x7);
0374         wire = ((iden >> 3) & 0x7f);
0375 
0376         meWheelOccup->Fill((float)wheel);
0377         meStationOccup->Fill((float)station);
0378         meSectorOccup->Fill((float)sector);
0379         meSuperLOccup->Fill((float)superlayer);
0380         meLayerOccup->Fill((float)layer);
0381         meWireOccup->Fill((float)wire);
0382 
0383         // Define a quantity to take into account station, splayer and layer being hit.
0384         path = (station - 1) * 40 + superlayer * 10 + layer;
0385         mePathMuon->Fill((float)path);
0386 
0387         // Define a quantity to take into chamber being hit.
0388         pathchamber = (wheel + 2) * 50 + (station - 1) * 12 + sector;
0389         meChamberOccup->Fill((float)pathchamber);
0390 
0391         /// Muon Momentum at MB1
0392         if (station == 1) {
0393           if (touch1 == 0) {
0394             mom1 = itHit->pabs();
0395             meMomentumMB1->Fill(mom1);
0396             touch1 = 1;
0397           }
0398         }
0399 
0400         /// Muon Momentum at MB4 & Loss of Muon Momentum in Iron (between MB1 and MB4)
0401         if (station == 4) {
0402           if (touch4 == 0) {
0403             mom4 = itHit->pabs();
0404             touch4 = 1;
0405             meMomentumMB4->Fill(mom4);
0406             if (touch1 == 1) {
0407               meLossMomIron->Fill(mom1 - mom4);
0408             }
0409           }
0410         }
0411 
0412         /// X-Local Coordinate vs Z-Local Coordinate
0413         meLocalXvsZ->Fill(itHit->localPosition().x(), itHit->localPosition().z());
0414 
0415         /// X-Local Coordinate vs Y-Local Coordinate
0416         meLocalXvsY->Fill(itHit->localPosition().x(), itHit->localPosition().y());
0417 
0418         /// Global Coordinates
0419 
0420         globposz = bsurf.toGlobal(itHit->localPosition()).z();
0421         globposeta = bsurf.toGlobal(itHit->localPosition()).eta();
0422         globposphi = bsurf.toGlobal(itHit->localPosition()).phi();
0423 
0424         radius = globposz * (1. + exp(-2. * globposeta)) / (1. - exp(-2. * globposeta));
0425 
0426         costeta = (1. - exp(-2. * globposeta)) / (1. + exp(-2. * globposeta));
0427         sinteta = 2. * exp(-globposeta) / (1. + exp(-2. * globposeta));
0428 
0429         /// Z-Global Coordinate vs X-Global Coordinate
0430         /// Y-Global Coordinate vs X-Global Coordinate
0431         globposx = radius * sinteta * cos(globposphi);
0432         globposy = radius * sinteta * sin(globposphi);
0433 
0434         meGlobalXvsZ->Fill(globposz, globposx);
0435         meGlobalXvsY->Fill(globposx, globposy);
0436 
0437         //  New Histos
0438         if (wheel == -2) {
0439           meGlobalXvsZWm2->Fill(globposz, globposx);
0440           meGlobalXvsYWm2->Fill(globposx, globposy);
0441         }
0442         if (wheel == -1) {
0443           meGlobalXvsZWm1->Fill(globposz, globposx);
0444           meGlobalXvsYWm1->Fill(globposx, globposy);
0445         }
0446         if (wheel == 0) {
0447           meGlobalXvsZW0->Fill(globposz, globposx);
0448           meGlobalXvsYW0->Fill(globposx, globposy);
0449         }
0450         if (wheel == 1) {
0451           meGlobalXvsZWp1->Fill(globposz, globposx);
0452           meGlobalXvsYWp1->Fill(globposx, globposy);
0453         }
0454         if (wheel == 2) {
0455           meGlobalXvsZWp2->Fill(globposz, globposx);
0456           meGlobalXvsYWp2->Fill(globposx, globposy);
0457         }
0458         //
0459         meHitRadius->Fill(radius);
0460         meCosTheta->Fill(costeta);
0461         meGlobalEta->Fill(globposeta);
0462         meGlobalPhi->Fill(globposphi);
0463       }
0464     } else {
0465       edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT")
0466           << "MuonDT PSimHit " << i << " is expected to be (det,subdet) = (" << dMuon << "," << sdMuonDT
0467           << "); value returned is: (" << detector << "," << subdetector << ")";
0468       continue;
0469     }
0470   }
0471 
0472   if (verbosity > 1) {
0473     eventout += "\n          Number of DT muon Hits collected:......... ";
0474     eventout += j;
0475   }
0476   meMuDTHits->Fill((float)nummu_DT);
0477 
0478   if (verbosity > 0)
0479     edm::LogInfo("MuonSimHitsValidAnalyzer::fillDT") << eventout << "\n";
0480   return;
0481 }