Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 /** derived from DTGeometryAnalyzer by Nicola Amapane
0002  *
0003  *  \author M. Maggi - INFN Bari
0004  */
0005 
0006 #include <memory>
0007 #include <fstream>
0008 #include <FWCore/Framework/interface/Frameworkfwd.h>
0009 
0010 #include <FWCore/Framework/interface/one/EDAnalyzer.h>
0011 #include <FWCore/Framework/interface/Event.h>
0012 #include <FWCore/Framework/interface/EventSetup.h>
0013 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0014 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0015 
0016 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
0017 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0018 #include <Geometry/CommonTopologies/interface/RectangularStripTopology.h>
0019 #include <Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h>
0020 
0021 #include <string>
0022 #include <cmath>
0023 #include <vector>
0024 #include <iomanip>
0025 #include <set>
0026 
0027 class RPCGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0028 public:
0029   RPCGeometryAnalyzer(const edm::ParameterSet& pset);
0030 
0031   ~RPCGeometryAnalyzer() override;
0032 
0033   void beginJob() override {}
0034   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0035   void endJob() override {}
0036 
0037   const std::string& myName() { return myName_; }
0038 
0039 private:
0040   const edm::ESGetToken<RPCGeometry, MuonGeometryRecord> tokRPC_;
0041   const int dashedLineWidth_;
0042   const std::string dashedLine_;
0043   const std::string myName_;
0044   std::ofstream ofos;
0045 };
0046 
0047 RPCGeometryAnalyzer::RPCGeometryAnalyzer(const edm::ParameterSet& /*iConfig*/)
0048     : tokRPC_{esConsumes<RPCGeometry, MuonGeometryRecord>(edm::ESInputTag{})},
0049       dashedLineWidth_(104),
0050       dashedLine_(std::string(dashedLineWidth_, '-')),
0051       myName_("RPCGeometryAnalyzer") {
0052   ofos.open("MytestOutput.out");
0053   edm::LogVerbatim("RPCGeometry") << "======================== Opening output file";
0054 }
0055 
0056 RPCGeometryAnalyzer::~RPCGeometryAnalyzer() {
0057   ofos.close();
0058   edm::LogVerbatim("RPCGeometry") << "======================== Closing output file";
0059 }
0060 
0061 void RPCGeometryAnalyzer::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0062   const RPCGeometry* pDD = &iSetup.getData(tokRPC_);
0063 
0064   ofos << myName() << ": Analyzer..." << std::endl;
0065   ofos << "start " << dashedLine_ << std::endl;
0066 
0067   ofos << " Geometry node for RPCGeom is  " << &(*pDD) << std::endl;
0068   ofos << " I have " << pDD->detTypes().size() << " detTypes" << std::endl;
0069   ofos << " I have " << pDD->detUnits().size() << " detUnits" << std::endl;
0070   ofos << " I have " << pDD->dets().size() << " dets" << std::endl;
0071   ofos << " I have " << pDD->rolls().size() << " rolls" << std::endl;
0072   ofos << " I have " << pDD->chambers().size() << " chambers" << std::endl;
0073 
0074   ofos << myName() << ": Begin iteration over geometry..." << std::endl;
0075   ofos << "iter " << dashedLine_ << std::endl;
0076 
0077   ofos << "\n  #     id(hex)      id(dec)                   "
0078           "  g(x=0)   g(y=0)   g(z=0)  g(z=-1)  g(z=+1)  Ns "
0079           "  phi(0)  phi(s1)  phi(sN)    dphi    dphi'      ds     off"
0080           "       uR       uL       lR       lL"
0081        << std::endl;
0082 
0083   int iRPCCHcount = 0;
0084 
0085   //  const double dPi = 3.14159265358;
0086   //   const double radToDeg = 180. / dPi; //@@ Where to get pi from?
0087 
0088   std::set<RPCDetId> sids;
0089   std::vector<LocalPoint> vlp;
0090   vlp.emplace_back(LocalPoint(-1, 0, 0));
0091   vlp.emplace_back(LocalPoint(0, 0, 0));
0092   vlp.emplace_back(LocalPoint(1, 0, 0));
0093   vlp.emplace_back(LocalPoint(0, -1, 0));
0094   vlp.emplace_back(LocalPoint(0, 0, 0));
0095   vlp.emplace_back(LocalPoint(0, 1, 0));
0096   vlp.emplace_back(LocalPoint(0, 0, -1));
0097   vlp.emplace_back(LocalPoint(0, 0, 0));
0098   vlp.emplace_back(LocalPoint(0, 0, 1));
0099 
0100   for (auto it : pDD->dets()) {
0101     //      //----------------------- RPCCHAMBER TEST -------------------------------------------------------
0102 
0103     if (dynamic_cast<const RPCChamber*>(it) != nullptr) {
0104       ++iRPCCHcount;
0105       const RPCChamber* ch = dynamic_cast<const RPCChamber*>(it);
0106 
0107       RPCDetId detId = ch->id();
0108       int idRaf = detId.rawId();
0109       //       const RPCRoll* rollRaf = ch->roll(1);
0110       ofos << "Num = " << iRPCCHcount << "  "
0111            << "RPCDet = " << idRaf << "  Num Rolls =" << ch->nrolls();
0112       //       "  "<<"Roll 1 = "<<(rollRaf->id()).rawId()<<std::endl;
0113 
0114       std::vector<const RPCRoll*> rollsRaf = (ch->rolls());
0115       for (auto& r : rollsRaf) {
0116         if (r->id().region() == 0) {
0117           ofos << "RPCDetId = " << r->id().rawId() << std::endl;
0118           ofos << "Region = " << r->id().region() << "  Ring = " << r->id().ring()
0119                << "  Station = " << r->id().station() << "  Sector = " << r->id().sector()
0120                << "  Layer = " << r->id().layer() << "  Subsector = " << r->id().subsector()
0121                << "  Roll = " << r->id().roll() << std::endl;
0122         }
0123       }
0124     }
0125     //_______________________________________________________________________________________________
0126 
0127     //      if( dynamic_cast< RPCRoll* >( *it ) != 0 ){
0128     //        ++icount;
0129 
0130     //        RPCRoll* roll = dynamic_cast< RPCRoll*>( *it );
0131 
0132     //        RPCDetId detId=roll->id();
0133     //        int id = detId(); // or detId.rawId()
0134 
0135     //        const StripTopology* top_ = dynamic_cast<const StripTopology*>
0136     //   (&roll->topology());
0137 
0138     //        if (sids.find(detId) != sids.end())
0139     //          st1 << "VERYBAD ";
0140     //          st1 << "GeomDetUnit is of type " << detId.det() << " and raw id = " << id;
0141 
0142     //          int iw = std::cout.width(); // save current width
0143     //          int ip = std::cout.precision(); // save current precision
0144 
0145     //          edm::LogVerbatim("RPCGeometry") << st1.str() << "Parameters of roll# "
0146     //                                          << std::setw( 4 ) << icount << std::endl;
0147     //        st2 << std::setw(12) << id << std::dec << std::setw(12)
0148     //            << id << std::dec << std::setw(iw) << detId;
0149 
0150     //        const BoundSurface& bSurface = roll->surface();
0151 
0152     //        LocalPoint  lCentre( 0., 0., 0. );
0153     //        GlobalPoint gCentre = bSurface.toGlobal( lCentre );
0154 
0155     //        LocalPoint  lCentre1( 0., 0., -1.);
0156     //        GlobalPoint gCentre1 = bSurface.toGlobal( lCentre1 );
0157 
0158     //        LocalPoint  lCentre2( 0., 0., 1.);
0159     //        GlobalPoint gCentre2 = bSurface.toGlobal( lCentre2 );
0160 
0161     //        double gx  =  gCentre.x();
0162     //        double gy  =  gCentre.y();
0163     //        double gz  =  gCentre.z();
0164     //        double gz1 =  gCentre1.z();
0165     //        double gz2 =  gCentre2.z();
0166     //        if ( fabs( gx )  < 1.e-06 ) gx = 0.;
0167     //        if ( fabs( gy )  < 1.e-06 ) gy = 0.;
0168     //        if ( fabs( gz )  < 1.e-06 ) gz = 0.;
0169     //        if ( fabs( gz1 ) < 1.e-06 ) gz1 = 0.;
0170     //        if ( fabs( gz2 ) < 1.e-06 ) gz2 = 0.;
0171 
0172     //        int now = 9;
0173     //        int nop = 5;
0174     //        edm::LogVerbatim("RPCGeometry") << st2.str() <<
0175     //       std::setw( now ) << std::setprecision( nop ) << gx <<
0176     //       std::setw( now ) << std::setprecision( nop ) << gy <<
0177     //       std::setw( now ) << std::setprecision( nop ) << gz <<
0178     //       std::setw( now ) << std::setprecision( nop ) << gz1 <<
0179     //       std::setw( now ) << std::setprecision( nop ) << gz2 << std::endl;
0180 
0181     //        // Global Phi of centre of RPCRoll
0182 
0183     //        double cphi = gCentre.phi();
0184     //        double cphiDeg = cphi * radToDeg;
0185     //        if( cphiDeg < 0. ) cphiDeg = cphiDeg + 360.;
0186 
0187     //        if ( fabs(cphiDeg) < 1.e-06 ) cphiDeg = 0.;
0188     //        //        int iphiDeg = static_cast<int>( cphiDeg );
0189     //        //    edm::LogVerbatim("RPCGeometry") << "phi(0,0,0) = " << iphiDeg << " degrees";
0190 
0191     //        int nStrips = roll->nstrips();
0192 
0193     //        if ( (detId.region()!=0 && detId.sector() == 1
0194     //       && detId.subsector() == 1
0195     //       && detId.roll() == 1
0196     //       && detId.station() == 1
0197     //       //&& detId.ring()==3
0198     //       && roll->type().name() == "REA1")
0199     //      ||
0200     //      (detId.region() == 0 && detId.sector() == 1
0201     //       && detId.station() == 1 && detId.layer() == 1
0202     //       && detId.roll() == 1
0203     //       && detId.ring() == 0)
0204     //      ) {
0205     //   edm::LogVerbatim("RPCGeometry") << "======================== Writing output file";
0206     //   ofos << "Forward Detector " << roll->type().name() << " " << detId.region() << " z :" << detId << std::endl;
0207     //   for (unsigned int i=0;i<vlp.size();i++){
0208     //     ofos << "lp=" << vlp[i] << " gp="<<roll->toGlobal(vlp[i]) << " pitch=" << roll->localPitch(vlp[i]);
0209     //     if ( (i+1)%3 == 0 ) {
0210     //       ofos << " " << std::endl;
0211     //     }
0212     //   }
0213     //   ofos << "            Navigating " << std::endl;
0214     //   LocalPoint edge1 = top_->localPosition(0.);
0215     //   LocalPoint edge2 = top_->localPosition((float)nStrips);
0216     //   float lsl1 = top_->localStripLength(edge1);
0217     //   float lsl2 = top_->localStripLength(edge2);
0218     //   ofos <<" Local Point edge1 = " << edge1 << " StripLength=" << lsl1 << std::endl;
0219     //   ofos <<" Local Point edge1 = " << edge2 << " StripLength=" << lsl2 << std::endl;
0220     //   float cslength = top_->localStripLength(lCentre);
0221     //   LocalPoint s1(0.,-cslength/2.,0.);
0222     //   LocalPoint s2(0.,cslength/2.,0.);
0223     //   float base1=top_->localPitch(s1)*nStrips;
0224     //   float base2=top_->localPitch(s2)*nStrips;
0225     //   //  LocalPoint  s11(-base1/2., -cslength/2,0.);
0226     //   //  LocalPoint  s12(base1/2., -cslength/2,0.);
0227     //   //  LocalPoint  s21(-base2/2., cslength/2,0.);
0228     //   //  LocalPoint  s22(base2/2.,  cslength/2,0.);
0229     //   ofos<<  "  First Base = "<<base1<<" Second Base ="<<base2<<std::endl;
0230     //        }
0231     // //        std::cout << "\nStrips =  "<<std::setw( 4 ) << nStrips<<"\n";
0232     //        for(int is=0;is<nStrips;is++){
0233     // //    std::cout <<"s="<<std::setw(3)<<is+1<<" pos="
0234     // //          <<roll->centreOfStrip(is+1);
0235     //   if ((is+1)%5==0){
0236     //     float str=is;
0237     // //      std::cout <<"s="<<std::setw(6)<<str<<" pos="
0238     // //            <<roll->centreOfStrip(str)<<" gpos="<<
0239     // //        roll->toGlobal(roll->centreOfStrip(str));
0240     // //      std::cout <<std::endl;
0241     //   }
0242   }
0243   ofos << std::endl;
0244 
0245   // //       double cstrip1  = layer->centerOfStrip(1).phi();
0246   // //       double cstripN  = layer->centerOfStrip(nStrips).phi();
0247 
0248   // //       double phiwid   = layer->geometry()->stripPhiPitch();
0249   // //       double stripwid = layer->geometry()->stripPitch();
0250   // //       double stripoff = layer->geometry()->stripOffset();
0251   // //       double phidif   = fabs(cstrip1 - cstripN);
0252 
0253   // //       // May have one strip at 180-epsilon and other at -180+epsilon
0254   // //       // If so the raw difference is 360-(phi extent of chamber)
0255   // //       // Want to reset that to (phi extent of chamber):
0256   // //       if ( phidif > dPi ) phidif = fabs(phidif - 2.*dPi);
0257   // //       double phiwid_check = phidif/double(nStrips-1);
0258 
0259   // //       // Clean up some stupid floating decimal aesthetics
0260   // //       cstrip1 = cstrip1 * radToDeg;
0261   // //       if ( cstrip1 < 0. ) cstrip1 = cstrip1 + 360.;
0262   // //       if ( fabs( cstrip1 ) < 1.e-06 ) cstrip1 = 0.;
0263   // //       cstripN = cstripN * radToDeg;
0264   // //       if ( cstripN < 0. ) cstripN = cstrip1 + 360.;
0265   // //       if ( fabs( cstripN ) < 1.e-06 ) cstripN = 0.;
0266 
0267   // //       if ( fabs( stripoff ) < 1.e-06 ) stripoff = 0.;
0268 
0269   // //       now = 9;
0270   // //       nop = 4;
0271   // //       ofos
0272   // //     << std::setw( now ) << std::setprecision( nop ) << cphiDeg
0273   // //     << std::setw( now ) << std::setprecision( nop ) << cstrip1
0274   // //     << std::setw( now ) << std::setprecision( nop ) << cstripN
0275   // //     << std::setw( now ) << std::setprecision( nop ) << phiwid
0276   // //     << std::setw( now ) << std::setprecision( nop ) << phiwid_check
0277   // //     << std::setw( now ) << std::setprecision( nop ) << stripwid
0278   // //     << std::setw( now ) << std::setprecision( nop ) << stripoff ;
0279   // //       //          << std::setw(8) << (layer->getOwner())->sector() ; //@@ No sector yet!
0280 
0281   // //       // Layer geometry:  layer corner phi's...
0282 
0283   // //       std::vector<float> parameters = layer->geometry()->parameters();
0284   // //       // these parameters are half-lengths, due to GEANT
0285   // //       float bottomEdge = parameters[0];
0286   // //       float topEdge    = parameters[1];
0287   // //       float thickness  = parameters[2];
0288   // //       float apothem    = parameters[3];
0289 
0290   // //       // first the face nearest the interaction
0291   // //       // get the other face by using positive thickness
0292   // //       LocalPoint upperRightLocal(topEdge, apothem, -thickness);
0293   // //       LocalPoint upperLeftLocal(-topEdge, apothem, -thickness);
0294   // //       LocalPoint lowerRightLocal(bottomEdge, -apothem, -thickness);
0295   // //       LocalPoint lowerLeftLocal(-bottomEdge, -apothem, -thickness);
0296 
0297   // //       GlobalPoint upperRightGlobal = bSurface.toGlobal(upperRightLocal);
0298   // //       GlobalPoint upperLeftGlobal  = bSurface.toGlobal(upperLeftLocal);
0299   // //       GlobalPoint lowerRightGlobal = bSurface.toGlobal(lowerRightLocal);
0300   // //       GlobalPoint lowerLeftGlobal  = bSurface.toGlobal(lowerLeftLocal);
0301 
0302   // //       float uRGp = radToDeg * upperRightGlobal.phi();
0303   // //       float uLGp = radToDeg * upperLeftGlobal.phi();
0304   // //       float lRGp = radToDeg * lowerRightGlobal.phi();
0305   // //       float lLGp = radToDeg * lowerLeftGlobal.phi();
0306 
0307   // //       now = 9;
0308   // //       ofos
0309   // //     << std::setw( now ) << uRGp
0310   // //     << std::setw( now ) << uLGp
0311   // //     << std::setw( now ) << lRGp
0312   // //     << std::setw( now ) << lLGp
0313   // //     << std::endl;
0314 
0315   // //       // Reset the values we changed
0316   // //       ofos << std::setprecision( ip ) << std::setw( iw );
0317 
0318   //        sids.insert(detId);
0319   //      }
0320   //    }
0321 
0322   // //   //   for (TrackingGeometry::DetTypeContainer::const_iterator it = pDD->detTypes().begin(); it != pDD->detTypes().end(); it ++){
0323   // //   //   }
0324 
0325   edm::LogVerbatim("RPCGeometry") << dashedLine_ << " end";
0326 }
0327 
0328 //define this as a plug-in
0329 #include <FWCore/Framework/interface/MakerMacros.h>
0330 DEFINE_FWK_MODULE(RPCGeometryAnalyzer);