EcalTestBeam/plugins/EcalTBMCInfoProducer.cc

0001 /*
0002  * \file EcalTBMCInfoProducer.cc
0003  *
0004  *
0005  */
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/Math/interface/Point3D.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/PluginManager/interface/ModuleDef.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0020
0021 #include "Geometry/EcalTestBeam/interface/EcalTBCrystalMap.h"
0022 #include "SimDataFormats/EcalTestBeam/interface/PEcalTBInfo.h"
0023 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0024
0025 #include "Math/GenVector/Rotation3D.h"
0026 #include <CLHEP/Random/RandFlat.h>
0027
0028 #include <fstream>
0029 #include <iostream>
0030 #include <memory>
0031 #include <vector>
0032
0033 class EcalTBMCInfoProducer : public edm::stream::EDProducer<> {
0034 public:
0035   /// Constructor
0036   explicit EcalTBMCInfoProducer(const edm::ParameterSet &ps);
0037
0038   /// Destructor
0039   ~EcalTBMCInfoProducer() override = default;
0040
0041   /// Produce digis out of raw data
0042   void produce(edm::Event &event, const edm::EventSetup &eventSetup) override;
0043
0044 private:
0045   const double fMinEta;
0046   const double fMaxEta;
0047   const double fMinPhi;
0048   const double fMaxPhi;
0049   const double beamEta;
0050   const double beamPhi;
0051   const double beamTheta;
0052   const double beamXoff;
0053   const double beamYoff;
0054
0055   const edm::EDGetTokenT<edm::HepMCProduct> GenVtxToken;
0056
0057   int crysNumber;
0058
0059   double partXhodo;
0060   double partYhodo;
0061
0062   std::unique_ptr<EcalTBCrystalMap> theTestMap;
0063
0064   std::unique_ptr<ROOT::Math::Rotation3D> fromCMStoTB;
0065 };
0066
0067 EcalTBMCInfoProducer::EcalTBMCInfoProducer(const edm::ParameterSet &ps)
0068     : fMinEta(ps.getParameter<double>("MinEta")),
0069       fMaxEta(ps.getParameter<double>("MaxEta")),
0070       fMinPhi(ps.getParameter<double>("MinPhi")),
0071       fMaxPhi(ps.getParameter<double>("MaxPhi")),
0072       beamEta((fMaxEta + fMinEta) / 2.),
0073       beamPhi((fMaxPhi + fMinPhi) / 2.),
0074       beamTheta(2.0 * atan(exp(-beamEta))),
0075       beamXoff(ps.getParameter<double>("BeamMeanX")),
0076       beamYoff(ps.getParameter<double>("BeamMeanX")),
0077       GenVtxToken(consumes<edm::HepMCProduct>(edm::InputTag("moduleLabelVtx", "source"))) {
0078   produces<PEcalTBInfo>();
0079
0080   edm::FileInPath CrystalMapFile = ps.getParameter<edm::FileInPath>("CrystalMapFile");
0081
0082   std::string fullMapName = CrystalMapFile.fullPath();
0083   theTestMap = std::make_unique<EcalTBCrystalMap>(fullMapName);
0084   crysNumber = 0;
0085
0086   double deltaEta = 999.;
0087   double deltaPhi = 999.;
0088   for (int cryIndex = 1; cryIndex <= EcalTBCrystalMap::NCRYSTAL; ++cryIndex) {
0089     double eta = 0;
0090     double phi = 0.;
0091     theTestMap->findCrystalAngles(cryIndex, eta, phi);
0092     if (fabs(beamEta - eta) < deltaEta && fabs(beamPhi - phi) < deltaPhi) {
0093       deltaEta = fabs(beamEta - eta);
0094       deltaPhi = fabs(beamPhi - phi);
0095       crysNumber = cryIndex;
0096     } else if (fabs(beamEta - eta) < deltaEta && fabs(beamPhi - phi) > deltaPhi) {
0097       if (fabs(beamPhi - phi) < 0.017) {
0098         deltaEta = fabs(beamEta - eta);
0099         deltaPhi = fabs(beamPhi - phi);
0100         crysNumber = cryIndex;
0101       }
0102     } else if (fabs(beamEta - eta) > deltaEta && fabs(beamPhi - phi) < deltaPhi) {
0103       if (fabs(beamEta - eta) < 0.017) {
0104         deltaEta = fabs(beamEta - eta);
0105         deltaPhi = fabs(beamPhi - phi);
0106         crysNumber = cryIndex;
0107       }
0108     }
0109   }
0110
0111   edm::LogVerbatim("EcalTBInfo") << "Initialize TB MC ECAL info producer with parameters: \n"
0112                                  << "Crystal map file:  " << CrystalMapFile << "\n"
0113                                  << "Beam average eta = " << beamEta << "\n"
0114                                  << "Beam average phi = " << beamPhi << "\n"
0115                                  << "Corresponding to crystal number = " << crysNumber << "\n"
0116                                  << "Beam X offset =    " << beamXoff << "\n"
0117                                  << "Beam Y offset =    " << beamYoff;
0118
0119   // rotation matrix to move from the CMS reference frame to the test beam one
0120
0121   double xx = -cos(beamTheta) * cos(beamPhi);
0122   double xy = -cos(beamTheta) * sin(beamPhi);
0123   double xz = sin(beamTheta);
0124
0125   double yx = sin(beamPhi);
0126   double yy = -cos(beamPhi);
0127   double yz = 0.;
0128
0129   double zx = sin(beamTheta) * cos(beamPhi);
0130   double zy = sin(beamTheta) * sin(beamPhi);
0131   double zz = cos(beamTheta);
0132
0133   fromCMStoTB = std::make_unique<ROOT::Math::Rotation3D>(xx, xy, xz, yx, yy, yz, zx, zy, zz);
0134
0135   // random number
0136   edm::Service<edm::RandomNumberGenerator> rng;
0137   if (!rng.isAvailable()) {
0138     throw cms::Exception("Configuration") << "EcalTBMCInfoProducer requires the RandomNumberGeneratorService\n"
0139                                              "which is not present in the configuration file.  You must add the "
0140                                              "service\n"
0141                                              "in the configuration file or remove the modules that require it.";
0142   }
0143 }
0144
0145 void EcalTBMCInfoProducer::produce(edm::Event &event, const edm::EventSetup &eventSetup) {
0146   edm::Service<edm::RandomNumberGenerator> rng;
0147   CLHEP::HepRandomEngine *engine = &rng->getEngine(event.streamID());
0148
0149   std::unique_ptr<PEcalTBInfo> product(new PEcalTBInfo());
0150
0151   // Fill the run information
0152
0153   product->setCrystal(crysNumber);
0154
0155   product->setBeamDirection(beamEta, beamPhi);
0156   product->setBeamOffset(beamXoff, beamYoff);
0157
0158   // Compute the event x,y vertex coordinates in the beam reference system
0159   // e.g. in the place orthogonal to the beam average direction
0160
0161   partXhodo = partYhodo = 0.;
0162
0163   const edm::Handle<edm::HepMCProduct> &GenEvt = event.getHandle(GenVtxToken);
0164
0165   const HepMC::GenEvent *Evt = GenEvt->GetEvent();
0166   HepMC::GenEvent::vertex_const_iterator Vtx = Evt->vertices_begin();
0167
0168   math::XYZPoint eventCMSVertex((*Vtx)->position().x(), (*Vtx)->position().y(), (*Vtx)->position().z());
0169
0170   LogDebug("EcalTBInfo") << "Generated vertex position = " << eventCMSVertex.x() << " " << eventCMSVertex.y() << " "
0171                          << eventCMSVertex.z();
0172
0173   math::XYZPoint eventTBVertex = (*fromCMStoTB) * eventCMSVertex;
0174
0175   LogDebug("EcalTBInfo") << "Rotated vertex position   = " << eventTBVertex.x() << " " << eventTBVertex.y() << " "
0176                          << eventTBVertex.z();
0177
0178   partXhodo = eventTBVertex.x();
0179   partYhodo = eventTBVertex.y();
0180
0181   product->setBeamPosition(partXhodo, partYhodo);
0182
0183   // Asynchronous phase shift
0184   double thisPhaseShift = CLHEP::RandFlat::shoot(engine);
0185
0186   product->setPhaseShift(thisPhaseShift);
0187   LogDebug("EcalTBInfo") << "Asynchronous Phaseshift = " << thisPhaseShift;
0188
0189   // store the object in the framework event
0190
0191   event.put(std::move(product));
0192 }
0193
0194 DEFINE_FWK_MODULE(EcalTBMCInfoProducer);