test/tools/MakeAngleLUT.cc

0001 //!
0002 //! This piece of code is obsolete and completely unused. It is kept as an
0003 //! example to show how to access certain geometry information.
0004 //!
0005
0006 #include <cassert>
0007 #include <cmath>
0008 #include <memory>
0009 #include <vector>
0010 #include <fstream>
0011 #include <iostream>
0012
0013 #include "TFile.h"
0014 #include "TTree.h"
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/Framework/interface/ConsumesCollector.h"
0024
0025 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0026 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0027 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0028 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0029 #include "MagneticField/Engine/interface/MagneticField.h"
0030 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0031
0032 //#include "DataFormats/MuonDetId/interface/CSCTriggerNumbering.h"
0033 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0034
0035 #include "L1Trigger/L1TMuon/interface/GeometryTranslator.h"
0036 #include "L1Trigger/L1TMuon/interface/MuonTriggerPrimitive.h"
0037 #include "L1Trigger/L1TMuon/interface/MuonTriggerPrimitiveFwd.h"
0038
0039 typedef L1TMuon::GeometryTranslator GeometryTranslator;
0040 typedef L1TMuon::TriggerPrimitive TriggerPrimitive;
0041 typedef L1TMuon::TriggerPrimitiveCollection TriggerPrimitiveCollection;
0042
0043 #include "helper.h"
0044
0045 class MakeAngleLUT : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0046 public:
0047   explicit MakeAngleLUT(const edm::ParameterSet&);
0048   ~MakeAngleLUT() override;
0049
0050 private:
0051   //virtual void beginJob();
0052   //virtual void endJob();
0053
0054   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0055   void endRun(const edm::Run&, const edm::EventSetup&) override;
0056
0057   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0058
0059   void generateLUTs();
0060
0061 private:
0062   GeometryTranslator geometry_translator_;
0063
0064   const edm::ParameterSet config_;
0065
0066   int verbose_;
0067
0068   std::string outfile_;
0069
0070   bool done_;
0071
0072   /// Event setup
0073 };
0074
0075 // _____________________________________________________________________________
0076 MakeAngleLUT::MakeAngleLUT(const edm::ParameterSet& iConfig)
0077     : geometry_translator_(consumesCollector()),
0078       config_(iConfig),
0079       verbose_(iConfig.getUntrackedParameter<int>("verbosity")),
0080       outfile_(iConfig.getParameter<std::string>("outfile")),
0081       done_(false) {
0082   assert(CSCConstants::KEY_CLCT_LAYER == CSCConstants::KEY_ALCT_LAYER);
0083 }
0084
0085 MakeAngleLUT::~MakeAngleLUT() {}
0086
0087 void MakeAngleLUT::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0088   geometry_translator_.checkAndUpdateGeometry(iSetup);
0089 }
0090
0091 void MakeAngleLUT::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0092
0093 void MakeAngleLUT::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0094   if (done_)
0095     return;
0096
0097   generateLUTs();
0098
0099   done_ = true;
0100   return;
0101 }
0102
0103 // _____________________________________________________________________________
0104 void MakeAngleLUT::generateLUTs() {
0105   const CSCGeometry& geocsc = geometry_translator_.getCSCGeometry();
0106   const RPCGeometry& georpc = geometry_translator_.getRPCGeometry();
0107   const GEMGeometry& geogem = geometry_translator_.getGEMGeometry();
0108   const MagneticField& magfield = geometry_translator_.getMagneticField();
0109
0110   auto average = [](double x, double y) { return 0.5 * (x + y); };
0111
0112   // Save z positions for ME1/1, ME1/2, ME1/3, ME2/2, ME3/2, ME4/2,
0113   //                      RE1/2, RE1/3, RE2/2, RE3/2, RE4/2,
0114   //                      GE1/1, GE1/2,
0115   std::vector<double> z_positions(13 * 2, 0.);
0116
0117   // ___________________________________________________________________________
0118   // CSC
0119
0120   for (const auto& it : geocsc.detUnits()) {
0121     const CSCLayer* layer = dynamic_cast<const CSCLayer*>(it);  // like GeomDetUnit
0122     assert(layer != nullptr);
0123     const CSCChamber* chamber = layer->chamber();  // like GeomDet
0124     assert(chamber != nullptr);
0125     const CSCDetId& cscDetId = chamber->id();
0126     //double zpos = chamber->surface().position().z();  // [cm]
0127     double zpos = chamber->layer(CSCConstants::KEY_ALCT_LAYER)->surface().position().z();  // [cm]
0128     //std::cout << "CSC: " << cscDetId.endcap() << " " << cscDetId.station() << " " << cscDetId.ring() << " " << cscDetId.chamber() << " " << cscDetId.layer() << " " << zpos << std::endl;
0129
0130     // Save the numbers
0131     if (cscDetId.endcap() == 1 && (cscDetId.chamber() == 1 || cscDetId.chamber() == 2)) {
0132       if (cscDetId.station() == 1 && cscDetId.ring() == 1) {
0133         if (cscDetId.chamber() == 2) {  // front
0134           z_positions[0] = zpos;
0135         } else if (cscDetId.chamber() == 1) {  // rear
0136           z_positions[1] = zpos;
0137         }
0138       } else if (cscDetId.station() == 1 && cscDetId.ring() == 2) {
0139         if (cscDetId.chamber() == 2) {  // front
0140           z_positions[2] = zpos;
0141         } else if (cscDetId.chamber() == 1) {  // rear
0142           z_positions[3] = zpos;
0143         }
0144       } else if (cscDetId.station() == 1 && cscDetId.ring() == 3) {
0145         if (cscDetId.chamber() == 2) {  // front
0146           z_positions[4] = zpos;
0147         } else if (cscDetId.chamber() == 1) {  // rear
0148           z_positions[5] = zpos;
0149         }
0150       } else if (cscDetId.station() == 2 && cscDetId.ring() == 2) {
0151         if (cscDetId.chamber() == 2) {  // front
0152           z_positions[6] = zpos;
0153         } else if (cscDetId.chamber() == 1) {  // rear
0154           z_positions[7] = zpos;
0155         }
0156       } else if (cscDetId.station() == 3 && cscDetId.ring() == 2) {
0157         if (cscDetId.chamber() == 1) {  // front
0158           z_positions[8] = zpos;
0159         } else if (cscDetId.chamber() == 2) {  // rear
0160           z_positions[9] = zpos;
0161         }
0162       } else if (cscDetId.station() == 4 && cscDetId.ring() == 2) {
0163         if (cscDetId.chamber() == 1) {  // front
0164           z_positions[10] = zpos;
0165         } else if (cscDetId.chamber() == 2) {  // rear
0166           z_positions[11] = zpos;
0167         }
0168       }
0169     }
0170   }  // end loop over CSC detUnits
0171
0172   // ___________________________________________________________________________
0173   // RPC
0174
0175   for (const auto& it : georpc.detUnits()) {
0176     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(it);  // like GeomDetUnit
0177     assert(roll != nullptr);
0178     //const RPCChamber* chamber = roll->chamber();  // like GeomDet
0179     //assert(chamber != nullptr);
0180     const RPCDetId& rpcDetId = roll->id();
0181     if (rpcDetId.region() == 0)  // skip barrel
0182       continue;
0183     //if (rpcDetId.region() == 0 || (rpcDetId.station() <= 2 && rpcDetId.ring() == 3))  // skip barrel, RE1/3, RE2/3
0184     //  continue;
0185     double zpos = roll->surface().position().z();  // [cm]
0186     //std::cout << "RPC: " << rpcDetId.region() << " " << rpcDetId.ring() << " " << rpcDetId.station() << " " << rpcDetId.sector() << " " << rpcDetId.layer() << " " << rpcDetId.subsector() << " " << rpcDetId.roll() << " " << zpos << std::endl;
0187
0188     // Save the numbers
0189     if (rpcDetId.region() == 1 && rpcDetId.sector() == 1 && rpcDetId.roll() == 1 &&
0190         (rpcDetId.subsector() == 1 || rpcDetId.subsector() == 2)) {
0191       if (rpcDetId.station() == 1 && rpcDetId.ring() == 2) {
0192         if (rpcDetId.subsector() == 2) {  // front
0193           z_positions[12] = zpos;
0194         } else if (rpcDetId.subsector() == 1) {  // rear
0195           z_positions[13] = zpos;
0196         }
0197       } else if (rpcDetId.station() == 1 && rpcDetId.ring() == 3) {
0198         if (rpcDetId.subsector() == 2) {  // front
0199           z_positions[14] = zpos;
0200         } else if (rpcDetId.subsector() == 1) {  // rear
0201           z_positions[15] = zpos;
0202         }
0203       } else if (rpcDetId.station() == 2 && rpcDetId.ring() == 2) {
0204         if (rpcDetId.subsector() == 2) {  // front
0205           z_positions[16] = zpos;
0206         } else if (rpcDetId.subsector() == 1) {  // rear
0207           z_positions[17] = zpos;
0208         }
0209       } else if (rpcDetId.station() == 3 && rpcDetId.ring() == 2) {
0210         if (rpcDetId.subsector() == 2) {  // front
0211           z_positions[18] = zpos;
0212         } else if (rpcDetId.subsector() == 1) {  // rear
0213           z_positions[19] = zpos;
0214         }
0215       } else if (rpcDetId.station() == 4 && rpcDetId.ring() == 2) {
0216         if (rpcDetId.subsector() == 2) {  // front
0217           z_positions[20] = zpos;
0218         } else if (rpcDetId.subsector() == 1) {  // rear
0219           z_positions[21] = zpos;
0220         }
0221       }
0222     }
0223   }  // end loop over RPC detUnits
0224
0225   // ___________________________________________________________________________
0226   // GEM
0227
0228   for (const auto& it : geogem.detUnits()) {
0229     const GEMEtaPartition* roll = dynamic_cast<const GEMEtaPartition*>(it);  // like GeomDetUnit
0230     assert(roll != nullptr);
0231     //const GEMChamber* chamber = roll->chamber();  // like GeomDet
0232     //assert(chamber != nullptr);
0233     const GEMDetId& gemDetId = roll->id();
0234     if (gemDetId.region() == 0)  // skip barrel
0235       continue;
0236     //double zpos = roll->surface().position().z();  // [cm]
0237     //std::cout << "GEM: " << gemDetId.region() << " " << gemDetId.ring() << " " << gemDetId.station() << " " << gemDetId.layer() << " " << gemDetId.chamber() << " " << gemDetId.roll() << " " << zpos << std::endl;
0238   }  // end loop over GEM detUnits
0239
0240   // ___________________________________________________________________________
0241   // Verbose
0242
0243   if (verbose_) {
0244     std::cout << "z positions:" << std::endl;
0245     for (const auto& zpos : z_positions) {
0246       std::cout << zpos << std::endl;
0247     }
0248     std::cout << std::endl;
0249
0250     std::cout << "common planes:" << std::endl;
0251     std::cout << average(z_positions[0], z_positions[1]) << std::endl;
0252     std::cout << average(z_positions[2], z_positions[3]) << std::endl;
0253     std::cout << average(z_positions[6], z_positions[7]) << std::endl;
0254     std::cout << average(z_positions[8], z_positions[9]) << std::endl;
0255     std::cout << average(z_positions[10], z_positions[11]) << std::endl;
0256     std::cout << std::endl;
0257   }
0258
0259   // Calculate the coefficients
0260   auto get_eta_bin_center = [](int b) {
0261     // nbinsx, xlow, xup = 2048, 1.1, 2.5
0262     double c = 1.1 + (2.5 - 1.1) / 2048. * (0.5 + static_cast<double>(b));
0263     return c;
0264   };
0265
0266   int num_z_bins = z_positions.size();
0267   int num_eta_bins = 2048;
0268
0269   std::vector<double> coefficients;
0270
0271   // Loop over z bins
0272   for (int iz = 0; iz < num_z_bins; ++iz) {
0273     // Find common plane
0274     double common_zpos = 0.;
0275     if (iz == 0 || iz == 1) {  // ME1/1
0276       common_zpos = average(z_positions[0], z_positions[1]);
0277     } else if (iz == 2 || iz == 3 || iz == 4 || iz == 5) {  // ME1/2, ME1/3
0278       common_zpos = average(z_positions[2], z_positions[3]);
0279     } else if (iz == 6 || iz == 7) {  // ME2/2
0280       common_zpos = average(z_positions[6], z_positions[7]);
0281     } else if (iz == 8 || iz == 9) {  // ME3/2
0282       common_zpos = average(z_positions[8], z_positions[9]);
0283     } else if (iz == 10 || iz == 11) {  // ME4/2
0284       common_zpos = average(z_positions[10], z_positions[11]);
0285     } else if (iz == 12 || iz == 13 || iz == 14 || iz == 15) {  // RE1/2, RE1/3
0286       common_zpos = average(z_positions[2], z_positions[3]);
0287     } else if (iz == 16 || iz == 17) {  // RE2/2
0288       common_zpos = average(z_positions[6], z_positions[7]);
0289     } else if (iz == 18 || iz == 19) {  // RE3/2
0290       common_zpos = average(z_positions[8], z_positions[9]);
0291     } else if (iz == 20 || iz == 21) {  // RE4/2
0292       common_zpos = average(z_positions[10], z_positions[11]);
0293     }
0294
0295     // Loop over eta bins
0296     for (int ieta = 0; ieta < num_eta_bins; ++ieta) {
0297       // Find magnetic field strength
0298       double zpos = z_positions.at(iz);
0299       double deltaZ = zpos - common_zpos;
0300       double eta = get_eta_bin_center(ieta);
0301       double cotTheta = std::sinh(eta);
0302       double r = zpos / cotTheta;
0303
0304       const GlobalPoint gp(r, 0, zpos);      // phi = 0
0305       double bz = magfield.inTesla(gp).z();  // [Tesla]
0306
0307       // Calculate the coefficient
0308       // coeff = deltaZ * (-0.5) * 0.003 * bz
0309       double coeff = deltaZ * (-0.5) * 0.003 * bz;
0310       coefficients.push_back(coeff);
0311     }
0312   }
0313   assert(coefficients.size() == (unsigned)(num_z_bins * num_eta_bins));
0314
0315   // Write
0316   {
0317     TFile* tfile = TFile::Open(outfile_.c_str(), "RECREATE");
0318     TTree* ttree = new TTree("tree", "tree");
0319     ttree->Branch("coefficients", &coefficients);
0320     ttree->Fill();
0321     tfile->Write();
0322     std::cout << "Wrote file: " << outfile_ << std::endl;
0323   }
0324 }
0325
0326 // DEFINE THIS AS A PLUG-IN
0327 #include "FWCore/Framework/interface/MakerMacros.h"
0328 DEFINE_FWK_MODULE(MakeAngleLUT);