Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:19

0001 // -*- C++ -*-
0002 //
0003 // Package:    RPCGEO
0004 // Class:      RPCGEO
0005 //
0006 /**\class RPCGEO RPCGEO.cc rpcgeo/RPCGEO/src/RPCGEO.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  pts/91
0015 //         Created:  Wed Sep 26 17:08:29 CEST 2007
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <fstream>
0022 #include <iomanip>
0023 
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 
0031 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0032 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0033 #include <DataFormats/RPCDigi/interface/RPCDigi.h>
0034 #include <DataFormats/RPCDigi/interface/RPCDigiCollection.h>
0035 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
0036 
0037 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0038 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
0039 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
0040 #include <Geometry/CommonTopologies/interface/RectangularStripTopology.h>
0041 #include <Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h>
0042 
0043 class RPCGEO : public edm::one::EDAnalyzer<> {
0044 public:
0045   explicit RPCGEO(const edm::ParameterSet&);
0046   ~RPCGEO() override;
0047 
0048   void beginJob() override {}
0049   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0050   void endJob() override {}
0051 
0052 private:
0053   const edm::ESGetToken<RPCGeometry, MuonGeometryRecord> tokRPC_;
0054   std::ofstream detidsstream;
0055   std::ofstream detailstream;
0056 };
0057 
0058 RPCGEO::RPCGEO(const edm::ParameterSet& /*iConfig*/)
0059     : tokRPC_{esConsumes<RPCGeometry, MuonGeometryRecord>(edm::ESInputTag{})} {
0060   detailstream.open("RPCGeometry.out");
0061   detidsstream.open("RPCDetIdLst.out");
0062   edm::LogVerbatim("RPCGeometry") << "Opening output files :: RPCGeometry.out and RPCDetIdLst.out";
0063 }
0064 
0065 RPCGEO::~RPCGEO() {
0066   detailstream.close();
0067   detidsstream.close();
0068   edm::LogVerbatim("RPCGeometry") << "Closing output files :: RPCGeometry.out and RPCDetIdLst.out";
0069 }
0070 
0071 // ------------ method called to for each event  ------------
0072 void RPCGEO::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0073   using namespace edm;
0074 
0075   detailstream << " Getting the RPC Geometry" << std::endl;
0076   const RPCGeometry* rpcGeo = &iSetup.getData(tokRPC_);
0077 
0078   int StripsInCMS = 0;
0079   int counterstripsBarrel = 0;
0080   int counterstripsEndCap = 0;
0081   int RollsInCMS = 0;
0082   int counterRollsBarrel = 0;
0083   int counterRollsMB1MB2MB3 = 0;
0084   int counterRollsMB4 = 0;
0085   int RB4Wm2 = 0;
0086   int RB4Wm1 = 0;
0087   int RB4W0 = 0;
0088   int RB4W1 = 0;
0089   int RB4W2 = 0;
0090   int counterRollsEndCap = 0;
0091   int ENDCAP[5][4];
0092   int ENDCAProll[5][4];
0093   int rollsNearDiskp3 = 0;
0094   int rollsNearDiskp2 = 0;
0095   float sumstripwbarrel = 0;
0096   float sumstripwendcap = 0;
0097   float areabarrel = 0;
0098   float areaendcap = 0;
0099 
0100   // Count chambers and rolls in the endcap
0101   // 4 Endcap disks with each 3 rings
0102   // count the chambers
0103   for (int i = 1; i < 5; i++) {
0104     for (int j = 1; j < 4; j++) {
0105       ENDCAP[i][j] = 0;
0106     }
0107   }
0108   // count the rolls
0109   for (int i = 1; i < 5; i++) {
0110     for (int j = 1; j < 4; j++) {
0111       ENDCAProll[i][j] = 0;
0112     }
0113   }
0114 
0115   for (TrackingGeometry::DetContainer::const_iterator it = rpcGeo->dets().begin(); it < rpcGeo->dets().end(); it++) {
0116     // The DetId can be a chamber or a roll
0117     // Consider only the chambers and ask the rolls belonging to this chamber later on
0118     // So this is a loop on the chambers
0119     if (dynamic_cast<const RPCChamber*>(*it) != nullptr) {
0120       const RPCChamber* ch = dynamic_cast<const RPCChamber*>(*it);
0121       std::vector<const RPCRoll*> rolls = (ch->rolls());
0122       //detailstream<<"RPC Chamber"<<ch->id()<<std::endl;
0123       if (ch->id().region() == 1) {
0124         switch (ch->id().station()) {
0125           case 1:
0126             switch (ch->id().ring()) {
0127               case 1:
0128                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0129                 break;
0130               case 2:
0131                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0132                 break;
0133               case 3:
0134                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0135                 break;
0136             }
0137             break;
0138           case 2:
0139             switch (ch->id().ring()) {
0140               case 1:
0141                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0142                 break;
0143               case 2:
0144                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0145                 break;
0146               case 3:
0147                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0148                 break;
0149             }
0150             break;
0151           case 3:
0152             switch (ch->id().ring()) {
0153               case 1:
0154                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0155                 break;
0156               case 2:
0157                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0158                 break;
0159               case 3:
0160                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0161                 break;
0162             }
0163             break;
0164           case 4:
0165             switch (ch->id().ring()) {
0166               case 1:
0167                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0168                 break;
0169               case 2:
0170                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0171                 break;
0172               case 3:
0173                 ++ENDCAP[ch->id().station()][ch->id().ring()];
0174                 break;
0175             }
0176             break;
0177         }
0178       }
0179       // ======================================
0180       // Loop over the rolls inside the chamber
0181       // ======================================
0182       for (std::vector<const RPCRoll*>::const_iterator r = rolls.begin(); r != rolls.end(); ++r) {
0183         RPCDetId rpcId = (*r)->id();
0184         int stripsinthisroll = (*r)->nstrips();
0185         RPCGeomServ rpcsrv(rpcId);
0186         ++RollsInCMS;
0187         // detailstream<<rpcId<<rpcsrv.name()<<" strips="<<stripsinthisroll<<std::endl;
0188         //detailstream<<rpcId<<" - "<<rpcsrv.name()<<" - "<<rpcsrv.shortname()<<std::endl;
0189         //detailstream<<rpcsrv.name()<<std::endl;
0190 
0191         // start RPC Barrel
0192         // ----------------
0193         if (rpcId.region() == 0) {
0194           //detailstream<<"Getting the RPC Topolgy"<<std::endl;
0195           const RectangularStripTopology* top_ = dynamic_cast<const RectangularStripTopology*>(&((*r)->topology()));
0196           float s1 = static_cast<float>(1) - 0.5;
0197           float sLast = static_cast<float>(stripsinthisroll) - 0.5;
0198           float stripl = top_->stripLength();
0199           float stripw = top_->pitch();
0200           areabarrel = areabarrel + stripl * stripw * stripsinthisroll;
0201           sumstripwbarrel = sumstripwbarrel + stripw * stripsinthisroll;
0202           std::string name = rpcsrv.name();
0203 
0204           // +++ Printout +++
0205           // ++++++++++++++++
0206           detidsstream << rpcId.rawId() << " \t " << std::setw(24) << name << " \t " << rpcId << std::endl;
0207           detailstream << "All Info " << rpcId.rawId() << " = " << std::setw(24) << name << " || ";
0208           detailstream << " strips: length =" << stripl << "[cm] & width =" << stripw
0209                        << "[cm] & number =" << stripsinthisroll;
0210           detailstream << " || area roll =" << stripl * stripw << "[cm^2] || area total barrel = " << areabarrel
0211                        << "[cm^2]" << std::endl;
0212           // ++++++++++++++++
0213 
0214           const BoundPlane& RPCSurface = (*r)->surface();
0215           GlobalPoint FirstStripCenterPointInGlobal = RPCSurface.toGlobal(top_->localPosition(s1));
0216           GlobalPoint LastStripCenterPointInGlobal = RPCSurface.toGlobal(top_->localPosition(sLast));
0217           double rpcphiFirst = FirstStripCenterPointInGlobal.barePhi();  //*180./3.141592;
0218           double rpcphiLast = LastStripCenterPointInGlobal.barePhi();    //*180./3.141592;
0219           //double rpcYFirst = FirstStripCenterPointInGlobal.y();
0220           //double rpcYLast  = LastStripCenterPointInGlobal.y();
0221           double diff = rpcphiLast - rpcphiFirst;
0222           double rollphi = (rpcphiFirst + rpcphiLast) * 0.5 * 180. / 3.141592;
0223           double orientation = diff / fabs(diff);
0224           int seg = rpcsrv.segment();
0225 
0226           // ++++++++++++++++
0227           // +++ Printout +++
0228           // ++++++++++++++++
0229           detailstream << std::setw(45) << name << " ||  phi=" << rollphi << " orientation=" << orientation
0230                        << " seg=" << seg << " phi first strip=" << rpcphiFirst << " phi last strip=" << rpcphiLast
0231                        << " || " << std::endl;
0232           detailstream << std::setw(45) << name << " ||  glob (X,Y,Z) first strip =("
0233                        << FirstStripCenterPointInGlobal.x() << ", " << FirstStripCenterPointInGlobal.y() << ","
0234                        << FirstStripCenterPointInGlobal.z() << ")[cm]"
0235                        << " glob (X,Y,Z) last strip =(" << LastStripCenterPointInGlobal.x() << ", "
0236                        << LastStripCenterPointInGlobal.y() << "," << LastStripCenterPointInGlobal.z() << ")[cm]"
0237                        << std::endl;
0238           // ++++++++++++++++
0239 
0240           ++counterRollsBarrel;
0241           if (rpcId.station() == 4) {
0242             ++counterRollsMB4;
0243             switch (rpcId.ring()) {
0244               case -2:
0245                 ++RB4Wm2;
0246                 break;
0247               case -1:
0248                 ++RB4Wm1;
0249                 break;
0250               case 0:
0251                 ++RB4W0;
0252                 break;
0253               case 1:
0254                 ++RB4W1;
0255                 break;
0256               case 2:
0257                 ++RB4W2;
0258                 break;
0259             }
0260           } else
0261             ++counterRollsMB1MB2MB3;
0262         }
0263         // end RPC Barrel loop
0264         // -------------------
0265 
0266         // start RPC Endcap
0267         // -------------------
0268         else {
0269           const TrapezoidalStripTopology* top_ = dynamic_cast<const TrapezoidalStripTopology*>(&((*r)->topology()));
0270           float s1 = static_cast<float>(1) - 0.5;
0271           float sLast = static_cast<float>(stripsinthisroll) - 0.5;
0272           float stripl = top_->stripLength();
0273           float stripw = top_->pitch();
0274           areaendcap = areaendcap + stripw * stripl * stripsinthisroll;
0275           sumstripwendcap = sumstripwendcap + stripw * stripsinthisroll;
0276           std::string name = rpcsrv.name();
0277 
0278           // +++ Printout +++
0279           // ++++++++++++++++
0280           detidsstream << rpcId.rawId() << " \t " << std::setw(24) << name << " \t " << rpcId << std::endl;
0281           detailstream << "All Info " << rpcId.rawId() << " = " << std::setw(24) << name << " || ";
0282           detailstream << " strips: length =" << stripl << "[cm] & width =" << stripw
0283                        << "[cm] & number =" << stripsinthisroll;
0284           detailstream << " || area roll =" << stripl * stripw << "[cm^2] || area total endcap = " << areaendcap
0285                        << "[cm^2]" << std::endl;
0286           // ++++++++++++++++
0287 
0288           const BoundPlane& RPCSurface = (*r)->surface();
0289           GlobalPoint FirstStripCenterPointInGlobal = RPCSurface.toGlobal(top_->localPosition(s1));
0290           GlobalPoint LastStripCenterPointInGlobal = RPCSurface.toGlobal(top_->localPosition(sLast));
0291           double rpcphiFirst = FirstStripCenterPointInGlobal.barePhi();  //*180./3.141592;
0292           double rpcphiLast = LastStripCenterPointInGlobal.barePhi();    //*180./3.141592;
0293           //double rpcYFirst = FirstStripCenterPointInGlobal.y();
0294           //double rpcYLast  = LastStripCenterPointInGlobal.y();
0295           double diff = rpcphiLast - rpcphiFirst;
0296           double rollphi = (rpcphiFirst + rpcphiLast) * 0.5 * 180. / 3.141592;
0297           double orientation = diff / fabs(diff);
0298           int seg = rpcsrv.segment();
0299           // ??? what does this mean ???
0300           if (seg == 19)
0301             orientation = orientation * -1;
0302           // ????
0303 
0304           // +++ Printout +++
0305           // ++++++++++++++++
0306           detailstream << std::setw(45) << name << " ||  phi=" << rollphi << " orientation=" << orientation
0307                        << " seg=" << seg << " phi first strip=" << rpcphiFirst << " phi last strip=" << rpcphiLast
0308                        << " || " << std::endl;
0309           detailstream << std::setw(45) << name << " ||  glob (X,Y,Z) first strip =("
0310                        << FirstStripCenterPointInGlobal.x() << ", " << FirstStripCenterPointInGlobal.y() << ","
0311                        << FirstStripCenterPointInGlobal.z() << ")[cm]"
0312                        << " glob (X,Y,Z) last strip =(" << LastStripCenterPointInGlobal.x() << ", "
0313                        << LastStripCenterPointInGlobal.y() << "," << LastStripCenterPointInGlobal.z() << ")[cm]";
0314           // ++++++++++++++++
0315           // cscphi = 2*3.1415926536+CenterPointCSCGlobal.barePhi():cscphi=CenterPointCSCGlobal.barePhi();
0316 
0317           // +++ Check the orientation +++
0318           // +++++++++++++++++++++++++++++
0319           bool orientation_ok = false;
0320           if (rpcId.station() == 1) {
0321             if (rpcId.ring() == 1 && orientation * rpcId.region() == 1.0) {
0322               orientation_ok = true;
0323             }
0324             if (rpcId.ring() == 2 && seg % 2 != 0 && orientation * rpcId.region() == 1.0) {
0325               orientation_ok = true;
0326             }
0327             if (rpcId.ring() == 2 && seg % 2 == 0 && orientation * rpcId.region() == -1.0) {
0328               orientation_ok = true;
0329             }
0330             if (rpcId.ring() == 3 && orientation * rpcId.region() == 1.0) {
0331               orientation_ok = true;
0332             }
0333           }
0334           if (rpcId.station() == 2 || rpcId.station() == 4) {
0335             if (orientation * rpcId.region() == -1.0) {
0336               orientation_ok = true;
0337             }
0338           }
0339           if (rpcId.station() == 3) {
0340             if (orientation * rpcId.region() == 1.0) {
0341               orientation_ok = true;
0342             }
0343           }
0344           if (orientation_ok)
0345             detailstream << " OK" << std::endl;
0346           else
0347             detailstream << " WRONG!!!" << std::endl;
0348           // +++++++++++++++++++++++++++++
0349 
0350           ++counterRollsEndCap;
0351           if (rpcId.region() == 1) {
0352             switch (rpcId.station()) {
0353               case 1:
0354                 switch (rpcId.ring()) {
0355                   case 1:
0356                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0357                     break;
0358                   case 2:
0359                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0360                     break;
0361                   case 3:
0362                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0363                     break;
0364                 }
0365                 break;
0366               case 2:
0367                 switch (rpcId.ring()) {
0368                   case 1:
0369                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0370                     break;
0371                   case 2:
0372                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0373                     break;
0374                   case 3:
0375                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0376                     break;
0377                 }
0378                 break;
0379               case 3:
0380                 switch (rpcId.ring()) {
0381                   case 1:
0382                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0383                     break;
0384                   case 2:
0385                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0386                     break;
0387                   case 3:
0388                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0389                     break;
0390                 }
0391                 break;
0392               case 4:
0393                 switch (rpcId.ring()) {
0394                   case 1:
0395                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0396                     break;
0397                   case 2:
0398                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0399                     break;
0400                   case 3:
0401                     ++ENDCAProll[rpcId.station()][rpcId.ring()];
0402                     break;
0403                 }
0404                 break;
0405             }
0406           }
0407         }
0408         // end RPC Endcap loop
0409         // -------------------
0410 
0411         // Particular Counter
0412         if (rpcId.region() == 1 && rpcId.station() == 2 &&
0413             (rpcId.sector() == 3 || rpcId.sector() == 4 || rpcId.sector() == 5)) {
0414           ++rollsNearDiskp2;
0415         }
0416         if (rpcId.region() == 1 && rpcId.station() == 3 &&
0417             (rpcId.sector() == 3 || rpcId.sector() == 4 || rpcId.sector() == 2)) {
0418           ++rollsNearDiskp3;
0419         }
0420         for (int strip = 1; strip <= stripsinthisroll; ++strip) {
0421           // LocalPoint lCentre=(*r)->centreOfStrip(strip);
0422           // const BoundSurface& bSurface = (*r)->surface();
0423           // GlobalPoint gCentre = bSurface.toGlobal(lCentre);
0424           // detailstream<<"Strip="<<strip<<" "<<gCentre.x()<<" "<<gCentre.y()<<" "<<gCentre.z()<<std::endl;
0425           ++StripsInCMS;
0426           if (rpcId.region() == 0)
0427             ++counterstripsBarrel;
0428           else
0429             ++counterstripsEndCap;
0430         }
0431       }
0432     }
0433   }
0434 
0435   detailstream << "Total Number of Strips in CMS=" << StripsInCMS << std::endl;
0436   detailstream << "Total Number of Rolls in CMS=" << RollsInCMS << std::endl;
0437 
0438   detailstream << "\n Total Number of Rolls in EndCap= " << counterRollsEndCap << std::endl;
0439   detailstream << "Total Number of Rolls in Barrel= " << counterRollsBarrel << std::endl;
0440 
0441   detailstream << "\n Total Number of Strips in Barrel= " << counterstripsBarrel << std::endl;
0442   detailstream << "Total Number of Strips in EndCap= " << counterstripsEndCap << std::endl;
0443 
0444   detailstream << "\n Total Number of Rolls in MB4= " << counterRollsMB4 << std::endl;
0445   detailstream << "Total Number of Rolls in MB1,MB2,MB3= " << counterRollsMB1MB2MB3 << std::endl;
0446 
0447   detailstream << "RB4 in the barrel:"
0448                << "\n Wheel -2 = " << RB4Wm2 << "\n Wheel -1 = " << RB4Wm1 << "\n Wheel 0 = " << RB4W0
0449                << "\n Wheel +1 = " << RB4W1 << "\n Wheel +2 = " << RB4W2 << std::endl;
0450 
0451   detailstream << "In the Endcaps we have the follow distribution of chambers:" << std::endl;
0452   for (int i = 1; i < 5; i++) {
0453     for (int j = 1; j < 4; j++) {
0454       detailstream << " Station " << i << " Ring " << j << " " << ENDCAP[i][j] << " Chambers" << std::endl;
0455     }
0456   }
0457 
0458   detailstream << "In the Endcaps we have the follow distribution of rolls:" << std::endl;
0459   for (int i = 1; i < 5; i++) {
0460     for (int j = 1; j < 4; j++) {
0461       detailstream << " Station " << i << " Ring " << j << " " << ENDCAProll[i][j] << " Roll" << std::endl;
0462     }
0463   }
0464 
0465   detailstream << "Rolls in Near Disk 2= " << rollsNearDiskp2 << std::endl;
0466   detailstream << "Rolls in Near Disk 3= " << rollsNearDiskp3 << std::endl;
0467 
0468   detailstream << "Average Strip Width in Barrel= " << sumstripwbarrel / counterstripsBarrel << std::endl;
0469   detailstream << "Average Strip Width in EndCap= " << sumstripwendcap / counterstripsEndCap << std::endl;
0470 
0471   detailstream << "Expected RMS Barrel= " << (sumstripwbarrel / counterstripsBarrel) / sqrt(12) << std::endl;
0472   detailstream << "Expected RMS EndCap= " << (sumstripwendcap / counterstripsEndCap) / sqrt(12) << std::endl;
0473 
0474   detailstream << "Area Covered in Barrel = " << areabarrel / 10000 << "m^2" << std::endl;
0475   detailstream << "Area Covered in EndCap = " << areaendcap / 10000 << "m^2" << std::endl;
0476 }
0477 
0478 //define this as a plug-in
0479 DEFINE_FWK_MODULE(RPCGEO);