Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:49 

0001 /** \class MuonDTDigis
0002  *  Analyse the the muon-drift-tubes digitizer.
0003  *
0004  *  \authors: R. Bellan
0005  */
0006 
0007 #include "MuonDTDigis.h"
0008 
0009 #include <iostream>
0010 #include <string>
0011 
0012 #include "TFile.h"
0013 
0014 #include "SimMuon/DTDigitizer/interface/Histograms.h"
0015 
0016 using namespace edm;
0017 using namespace std;
0018 
0019 MuonDTDigis::MuonDTDigis(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0020   // ----------------------
0021   // Get the debug parameter for verbose output
0022   verbose_ = pset.getUntrackedParameter<bool>("verbose", false);
0023 
0024   // the name of the Digi collection
0025   SimHitToken_ = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("SimHitLabel"));
0026 
0027   // the name of the Digi collection
0028   DigiToken_ = consumes<DTDigiCollection>(pset.getParameter<edm::InputTag>(("DigiLabel")));
0029 
0030   hDigis_global = std::make_unique<hDigis>("Global");
0031   hDigis_W0 = std::make_unique<hDigis>("Wheel0");
0032   hDigis_W1 = std::make_unique<hDigis>("Wheel1");
0033   hDigis_W2 = std::make_unique<hDigis>("Wheel2");
0034   hAllHits = std::make_unique<hHits>("AllHits");
0035 }
0036 
0037 MuonDTDigis::~MuonDTDigis() {
0038   if (verbose_)
0039     cout << "[MuonDTDigis] Destructor called" << endl;
0040 }
0041 
0042 void MuonDTDigis::bookHistograms(DQMStore::IBooker &iBooker,
0043                                  edm::Run const &iRun,
0044                                  edm::EventSetup const & /* iSetup */) {
0045   meDigiTimeBox_ = nullptr;
0046   meDigiTimeBox_wheel2m_ = nullptr;
0047   meDigiTimeBox_wheel1m_ = nullptr;
0048   meDigiTimeBox_wheel0_ = nullptr;
0049   meDigiTimeBox_wheel1p_ = nullptr;
0050   meDigiTimeBox_wheel2p_ = nullptr;
0051   meDigiEfficiency_ = nullptr;
0052   meDigiEfficiencyMu_ = nullptr;
0053   meDoubleDigi_ = nullptr;
0054   meSimvsDigi_ = nullptr;
0055   meWire_DoubleDigi_ = nullptr;
0056 
0057   meMB1_sim_occup_ = nullptr;
0058   meMB1_digi_occup_ = nullptr;
0059   meMB2_sim_occup_ = nullptr;
0060   meMB2_digi_occup_ = nullptr;
0061   meMB3_sim_occup_ = nullptr;
0062   meMB3_digi_occup_ = nullptr;
0063   meMB4_sim_occup_ = nullptr;
0064   meMB4_digi_occup_ = nullptr;
0065 
0066   meDigiHisto_ = nullptr;
0067 
0068   // ----------------------
0069   // We go
0070   // ----------------------
0071 
0072   Char_t histo_n[100];
0073   Char_t histo_t[100];
0074 
0075   iBooker.setCurrentFolder("MuonDTDigisV/DTDigiValidationTask");
0076 
0077   sprintf(histo_n, "DigiTimeBox");
0078   sprintf(histo_t, "Digi Time Box");
0079   meDigiTimeBox_ = iBooker.book1D(histo_n, histo_t, 1536, 0, 1200);
0080 
0081   sprintf(histo_n, "DigiTimeBox_wheel2m");
0082   sprintf(histo_t, "Digi Time Box wheel -2");
0083   meDigiTimeBox_wheel2m_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
0084 
0085   sprintf(histo_n, "DigiTimeBox_wheel1m");
0086   sprintf(histo_t, "Digi Time Box wheel -1");
0087   meDigiTimeBox_wheel1m_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
0088 
0089   sprintf(histo_n, "DigiTimeBox_wheel0");
0090   sprintf(histo_t, "Digi Time Box wheel 0");
0091   meDigiTimeBox_wheel0_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
0092 
0093   sprintf(histo_n, "DigiTimeBox_wheel1p");
0094   sprintf(histo_t, "Digi Time Box wheel 1");
0095   meDigiTimeBox_wheel1p_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
0096 
0097   sprintf(histo_n, "DigiTimeBox_wheel2p");
0098   sprintf(histo_t, "Digi Time Box wheel 2");
0099   meDigiTimeBox_wheel2p_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
0100 
0101   sprintf(histo_n, "DigiEfficiencyMu");
0102   sprintf(histo_t, "Ratio (#Digis Mu)/(#SimHits Mu)");
0103   meDigiEfficiencyMu_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);
0104 
0105   sprintf(histo_n, "DigiEfficiency");
0106   sprintf(histo_t, "Ratio (#Digis)/(#SimHits)");
0107   meDigiEfficiency_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);
0108 
0109   sprintf(histo_n, "Number_Digi_per_layer");
0110   sprintf(histo_t, "Number_Digi_per_layer");
0111   meDoubleDigi_ = iBooker.book1D(histo_n, histo_t, 10, 0., 10.);
0112 
0113   sprintf(histo_n, "Number_simhit_vs_digi");
0114   sprintf(histo_t, "Number_simhit_vs_digi");
0115   meSimvsDigi_ = iBooker.book2D(histo_n, histo_t, 100, 0., 140., 100, 0., 140.);
0116 
0117   sprintf(histo_n, "Wire_Number_with_double_Digi");
0118   sprintf(histo_t, "Wire_Number_with_double_Digi");
0119   meWire_DoubleDigi_ = iBooker.book1D(histo_n, histo_t, 100, 0., 100.);
0120 
0121   sprintf(histo_n, "Simhit_occupancy_MB1");
0122   sprintf(histo_t, "Simhit_occupancy_MB1");
0123   meMB1_sim_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55.);
0124 
0125   sprintf(histo_n, "Digi_occupancy_MB1");
0126   sprintf(histo_t, "Digi_occupancy_MB1");
0127   meMB1_digi_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55.);
0128 
0129   sprintf(histo_n, "Simhit_occupancy_MB2");
0130   sprintf(histo_t, "Simhit_occupancy_MB2");
0131   meMB2_sim_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63.);
0132 
0133   sprintf(histo_n, "Digi_occupancy_MB2");
0134   sprintf(histo_t, "Digi_occupancy_MB2");
0135   meMB2_digi_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63.);
0136 
0137   sprintf(histo_n, "Simhit_occupancy_MB3");
0138   sprintf(histo_t, "Simhit_occupancy_MB3");
0139   meMB3_sim_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75.);
0140 
0141   sprintf(histo_n, "Digi_occupancy_MB3");
0142   sprintf(histo_t, "Digi_occupancy_MB3");
0143   meMB3_digi_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75.);
0144 
0145   sprintf(histo_n, "Simhit_occupancy_MB4");
0146   sprintf(histo_t, "Simhit_occupancy_MB4");
0147   meMB4_sim_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99.);
0148 
0149   sprintf(histo_n, "Digi_occupancy_MB4");
0150   sprintf(histo_t, "Digi_occupancy_MB4");
0151   meMB4_digi_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99.);
0152 
0153   // Begona's option
0154   char stringcham[40];
0155   for (int slnum = 1; slnum < 62; ++slnum) {
0156     sprintf(stringcham, "DigiTimeBox_slid_%d", slnum);
0157     meDigiHisto_ = iBooker.book1D(stringcham, stringcham, 100, 0, 1200);
0158     meDigiTimeBox_SL_.push_back(meDigiHisto_);
0159   }
0160 }
0161 
0162 void MuonDTDigis::analyze(const Event &event, const EventSetup &eventSetup) {
0163   if (verbose_)
0164     cout << "--- [MuonDTDigis] Analysing Event: #Run: " << event.id().run() << " #Event: " << event.id().event()
0165          << endl;
0166 
0167   // Get the DT Geometry
0168   muonGeom = &eventSetup.getData(muonGeomToken_);
0169 
0170   // Get the Digi collection from the event
0171   Handle<DTDigiCollection> dtDigis;
0172   event.getByToken(DigiToken_, dtDigis);
0173 
0174   // Get simhits
0175   Handle<PSimHitContainer> simHits;
0176   event.getByToken(SimHitToken_, simHits);
0177 
0178   int num_mudigis;
0179   int num_musimhits;
0180   int num_digis;
0181   int num_digis_layer;
0182   int cham_num;
0183   int wire_touched;
0184   num_digis = 0;
0185   num_mudigis = 0;
0186   num_musimhits = 0;
0187   DTWireIdMap wireMap;
0188 
0189   for (vector<PSimHit>::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++) {
0190     // Create the id of the wire, the simHits in the DT known also the wireId
0191     DTWireId wireId(hit->detUnitId());
0192     //   cout << " PSimHits wire id " << wireId << " part type " <<
0193     //   hit->particleType() << endl;
0194 
0195     // Fill the map
0196     wireMap[wireId].push_back(&(*hit));
0197 
0198     LocalPoint entryP = hit->entryPoint();
0199     LocalPoint exitP = hit->exitPoint();
0200     int partType = hit->particleType();
0201     if (abs(partType) == 13)
0202       num_musimhits++;
0203 
0204     if (wireId.station() == 1 && abs(partType) == 13)
0205       meMB1_sim_occup_->Fill(wireId.wire());
0206     if (wireId.station() == 2 && abs(partType) == 13)
0207       meMB2_sim_occup_->Fill(wireId.wire());
0208     if (wireId.station() == 3 && abs(partType) == 13)
0209       meMB3_sim_occup_->Fill(wireId.wire());
0210     if (wireId.station() == 4 && abs(partType) == 13)
0211       meMB4_sim_occup_->Fill(wireId.wire());
0212 
0213     float path = (exitP - entryP).mag();
0214     float path_x = fabs((exitP - entryP).x());
0215 
0216     hAllHits->Fill(entryP.x(),
0217                    exitP.x(),
0218                    entryP.y(),
0219                    exitP.y(),
0220                    entryP.z(),
0221                    exitP.z(),
0222                    path,
0223                    path_x,
0224                    partType,
0225                    hit->processType(),
0226                    hit->pabs());
0227   }
0228 
0229   //  cout << "num muon simhits " << num_musimhits << endl;
0230 
0231   DTDigiCollection::DigiRangeIterator detUnitIt;
0232   for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
0233     const DTLayerId &id = (*detUnitIt).first;
0234     const DTDigiCollection::Range &range = (*detUnitIt).second;
0235 
0236     num_digis_layer = 0;
0237     cham_num = 0;
0238     wire_touched = 0;
0239 
0240     // Loop over the digis of this DetUnit
0241     for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0242       //   cout<<" Wire: "<<(*digiIt).wire()<<endl
0243       //  <<" digi time (ns): "<<(*digiIt).time()<<endl;
0244 
0245       num_digis++;
0246       num_digis_layer++;
0247       if (num_digis_layer > 1) {
0248         if ((*digiIt).wire() == wire_touched) {
0249           meWire_DoubleDigi_->Fill((*digiIt).wire());
0250           //      cout << "old & new wire " << wire_touched << " " <<
0251           //      (*digiIt).wire() << endl;
0252         }
0253       }
0254       wire_touched = (*digiIt).wire();
0255 
0256       meDigiTimeBox_->Fill((*digiIt).time());
0257       if (id.wheel() == -2)
0258         meDigiTimeBox_wheel2m_->Fill((*digiIt).time());
0259       if (id.wheel() == -1)
0260         meDigiTimeBox_wheel1m_->Fill((*digiIt).time());
0261       if (id.wheel() == 0)
0262         meDigiTimeBox_wheel0_->Fill((*digiIt).time());
0263       if (id.wheel() == 1)
0264         meDigiTimeBox_wheel1p_->Fill((*digiIt).time());
0265       if (id.wheel() == 2)
0266         meDigiTimeBox_wheel2p_->Fill((*digiIt).time());
0267 
0268       //   Superlayer number and fill histo with digi timebox
0269 
0270       cham_num = (id.wheel() + 2) * 12 + (id.station() - 1) * 3 + id.superlayer();
0271       //   cout << " Histo number " << cham_num << endl;
0272 
0273       meDigiTimeBox_SL_[cham_num]->Fill((*digiIt).time());
0274 
0275       //    cout << " size de digis " << (*digiIt).size() << endl;
0276 
0277       DTWireId wireId(id, (*digiIt).wire());
0278       if (wireId.station() == 1)
0279         meMB1_digi_occup_->Fill((*digiIt).wire());
0280       if (wireId.station() == 2)
0281         meMB2_digi_occup_->Fill((*digiIt).wire());
0282       if (wireId.station() == 3)
0283         meMB3_digi_occup_->Fill((*digiIt).wire());
0284       if (wireId.station() == 4)
0285         meMB4_digi_occup_->Fill((*digiIt).wire());
0286 
0287       int mu = 0;
0288       float theta = 0;
0289 
0290       for (vector<const PSimHit *>::iterator hit = wireMap[wireId].begin(); hit != wireMap[wireId].end(); hit++)
0291         if (abs((*hit)->particleType()) == 13) {
0292           theta = atan((*hit)->momentumAtEntry().x() / (-(*hit)->momentumAtEntry().z())) * 180 / M_PI;
0293           //      cout<<"momentum x: "<<(*hit)->momentumAtEntry().x()<<endl
0294           //          <<"momentum z: "<<(*hit)->momentumAtEntry().z()<<endl
0295           //          <<"atan: "<<theta<<endl;
0296           mu++;
0297         }
0298 
0299       if (mu)
0300         num_mudigis++;
0301 
0302       if (mu && theta) {
0303         hDigis_global->Fill((*digiIt).time(), theta, id.superlayer());
0304         // filling digi histos for wheel and for RZ and RPhi
0305         WheelHistos(id.wheel())->Fill((*digiIt).time(), theta, id.superlayer());
0306       }
0307 
0308     }  // for digis in layer
0309 
0310     meDoubleDigi_->Fill((float)num_digis_layer);
0311 
0312   }  // for layers
0313 
0314   // cout << "num_digis " << num_digis << "mu digis " << num_mudigis << endl;
0315 
0316   if (num_musimhits != 0) {
0317     meDigiEfficiencyMu_->Fill((float)num_mudigis / (float)num_musimhits);
0318     meDigiEfficiency_->Fill((float)num_digis / (float)num_musimhits);
0319   }
0320 
0321   meSimvsDigi_->Fill((float)num_musimhits, (float)num_digis);
0322   //  cout<<"--------------"<<endl;
0323 }
0324 
0325 hDigis *MuonDTDigis::WheelHistos(int wheel) {
0326   switch (abs(wheel)) {
0327     case 0:
0328       return hDigis_W0.get();
0329 
0330     case 1:
0331       return hDigis_W1.get();
0332 
0333     case 2:
0334       return hDigis_W2.get();
0335 
0336     default:
0337       return nullptr;
0338   }
0339 }
0340 #include "DQMServices/Core/interface/DQMStore.h"
0341 #include "FWCore/Framework/interface/MakerMacros.h"
0342 #include "FWCore/PluginManager/interface/ModuleDef.h"
0343 
0344 DEFINE_FWK_MODULE(MuonDTDigis);
[ 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!