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
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
0035 DTHitsToken_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("DTHitsSrc"));
0036
0037
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
0052
0053 << " DTHitsSrc = " << l.module << ":" << l.productInstance
0054 << "\n"
0055
0056
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& ) {
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
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
0257 ++count;
0258
0259
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
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
0310 edm::PSimHitContainer::const_iterator itHit;
0311
0312
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
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
0336 int i = 0, j = 0;
0337 for (itHit = MuonDTContainer->begin(); itHit != MuonDTContainer->end(); ++itHit) {
0338 ++i;
0339
0340
0341 DetId theDetUnitId(itHit->detUnitId());
0342 int detector = theDetUnitId.det();
0343 int subdetector = theDetUnitId.subdetId();
0344
0345
0346 if ((detector == dMuon) && (subdetector == sdMuonDT)) {
0347
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
0358 const BoundPlane& bsurf = theDet->surface();
0359
0360
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
0384 path = (station - 1) * 40 + superlayer * 10 + layer;
0385 mePathMuon->Fill((float)path);
0386
0387
0388 pathchamber = (wheel + 2) * 50 + (station - 1) * 12 + sector;
0389 meChamberOccup->Fill((float)pathchamber);
0390
0391
0392 if (station == 1) {
0393 if (touch1 == 0) {
0394 mom1 = itHit->pabs();
0395 meMomentumMB1->Fill(mom1);
0396 touch1 = 1;
0397 }
0398 }
0399
0400
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
0413 meLocalXvsZ->Fill(itHit->localPosition().x(), itHit->localPosition().z());
0414
0415
0416 meLocalXvsY->Fill(itHit->localPosition().x(), itHit->localPosition().y());
0417
0418
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
0430
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
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 }