Back to home page

Project CMSSW displayed by LXR

 
 

    


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