CMS 3D CMS Logo

EcalTBMCInfoProducer.cc
Go to the documentation of this file.
1 /*
2  * \file EcalTBMCInfoProducer.cc
3  *
4  *
5  */
6 
9 
20 
24 
25 #include "Math/GenVector/Rotation3D.h"
26 #include <CLHEP/Random/RandFlat.h>
27 
28 #include <fstream>
29 #include <iostream>
30 #include <memory>
31 #include <vector>
32 
34 public:
36  explicit EcalTBMCInfoProducer(const edm::ParameterSet &ps);
37 
39  ~EcalTBMCInfoProducer() override = default;
40 
42  void produce(edm::Event &event, const edm::EventSetup &eventSetup) override;
43 
44 private:
45  const double fMinEta;
46  const double fMaxEta;
47  const double fMinPhi;
48  const double fMaxPhi;
49  const double beamEta;
50  const double beamPhi;
51  const double beamTheta;
52  const double beamXoff;
53  const double beamYoff;
54 
56 
58 
59  double partXhodo;
60  double partYhodo;
61 
62  std::unique_ptr<EcalTBCrystalMap> theTestMap;
63 
64  std::unique_ptr<ROOT::Math::Rotation3D> fromCMStoTB;
65 };
66 
68  : fMinEta(ps.getParameter<double>("MinEta")),
69  fMaxEta(ps.getParameter<double>("MaxEta")),
70  fMinPhi(ps.getParameter<double>("MinPhi")),
71  fMaxPhi(ps.getParameter<double>("MaxPhi")),
72  beamEta((fMaxEta + fMinEta) / 2.),
73  beamPhi((fMaxPhi + fMinPhi) / 2.),
74  beamTheta(2.0 * atan(exp(-beamEta))),
75  beamXoff(ps.getParameter<double>("BeamMeanX")),
76  beamYoff(ps.getParameter<double>("BeamMeanX")),
77  GenVtxToken(consumes<edm::HepMCProduct>(edm::InputTag("moduleLabelVtx", "source"))) {
78  produces<PEcalTBInfo>();
79 
81 
82  std::string fullMapName = CrystalMapFile.fullPath();
83  theTestMap = std::make_unique<EcalTBCrystalMap>(fullMapName);
84  crysNumber = 0;
85 
86  double deltaEta = 999.;
87  double deltaPhi = 999.;
88  for (int cryIndex = 1; cryIndex <= EcalTBCrystalMap::NCRYSTAL; ++cryIndex) {
89  double eta = 0;
90  double phi = 0.;
91  theTestMap->findCrystalAngles(cryIndex, eta, phi);
92  if (fabs(beamEta - eta) < deltaEta && fabs(beamPhi - phi) < deltaPhi) {
93  deltaEta = fabs(beamEta - eta);
94  deltaPhi = fabs(beamPhi - phi);
95  crysNumber = cryIndex;
96  } else if (fabs(beamEta - eta) < deltaEta && fabs(beamPhi - phi) > deltaPhi) {
97  if (fabs(beamPhi - phi) < 0.017) {
98  deltaEta = fabs(beamEta - eta);
99  deltaPhi = fabs(beamPhi - phi);
100  crysNumber = cryIndex;
101  }
102  } else if (fabs(beamEta - eta) > deltaEta && fabs(beamPhi - phi) < deltaPhi) {
103  if (fabs(beamEta - eta) < 0.017) {
104  deltaEta = fabs(beamEta - eta);
105  deltaPhi = fabs(beamPhi - phi);
106  crysNumber = cryIndex;
107  }
108  }
109  }
110 
111  edm::LogVerbatim("EcalTBInfo") << "Initialize TB MC ECAL info producer with parameters: \n"
112  << "Crystal map file: " << CrystalMapFile << "\n"
113  << "Beam average eta = " << beamEta << "\n"
114  << "Beam average phi = " << beamPhi << "\n"
115  << "Corresponding to crystal number = " << crysNumber << "\n"
116  << "Beam X offset = " << beamXoff << "\n"
117  << "Beam Y offset = " << beamYoff;
118 
119  // rotation matrix to move from the CMS reference frame to the test beam one
120 
121  double xx = -cos(beamTheta) * cos(beamPhi);
122  double xy = -cos(beamTheta) * sin(beamPhi);
123  double xz = sin(beamTheta);
124 
125  double yx = sin(beamPhi);
126  double yy = -cos(beamPhi);
127  double yz = 0.;
128 
129  double zx = sin(beamTheta) * cos(beamPhi);
130  double zy = sin(beamTheta) * sin(beamPhi);
131  double zz = cos(beamTheta);
132 
133  fromCMStoTB = std::make_unique<ROOT::Math::Rotation3D>(xx, xy, xz, yx, yy, yz, zx, zy, zz);
134 
135  // random number
137  if (!rng.isAvailable()) {
138  throw cms::Exception("Configuration") << "EcalTBMCInfoProducer requires the RandomNumberGeneratorService\n"
139  "which is not present in the configuration file. You must add the "
140  "service\n"
141  "in the configuration file or remove the modules that require it.";
142  }
143 }
144 
147  CLHEP::HepRandomEngine *engine = &rng->getEngine(event.streamID());
148 
149  std::unique_ptr<PEcalTBInfo> product(new PEcalTBInfo());
150 
151  // Fill the run information
152 
153  product->setCrystal(crysNumber);
154 
155  product->setBeamDirection(beamEta, beamPhi);
156  product->setBeamOffset(beamXoff, beamYoff);
157 
158  // Compute the event x,y vertex coordinates in the beam reference system
159  // e.g. in the place orthogonal to the beam average direction
160 
161  partXhodo = partYhodo = 0.;
162 
163  const edm::Handle<edm::HepMCProduct> &GenEvt = event.getHandle(GenVtxToken);
164 
165  const HepMC::GenEvent *Evt = GenEvt->GetEvent();
166  HepMC::GenEvent::vertex_const_iterator Vtx = Evt->vertices_begin();
167 
168  math::XYZPoint eventCMSVertex((*Vtx)->position().x(), (*Vtx)->position().y(), (*Vtx)->position().z());
169 
170  LogDebug("EcalTBInfo") << "Generated vertex position = " << eventCMSVertex.x() << " " << eventCMSVertex.y() << " "
171  << eventCMSVertex.z();
172 
173  math::XYZPoint eventTBVertex = (*fromCMStoTB) * eventCMSVertex;
174 
175  LogDebug("EcalTBInfo") << "Rotated vertex position = " << eventTBVertex.x() << " " << eventTBVertex.y() << " "
176  << eventTBVertex.z();
177 
178  partXhodo = eventTBVertex.x();
179  partYhodo = eventTBVertex.y();
180 
181  product->setBeamPosition(partXhodo, partYhodo);
182 
183  // Asynchronous phase shift
184  double thisPhaseShift = CLHEP::RandFlat::shoot(engine);
185 
186  product->setPhaseShift(thisPhaseShift);
187  LogDebug("EcalTBInfo") << "Asynchronous Phaseshift = " << thisPhaseShift;
188 
189  // store the object in the framework event
190 
191  event.put(std::move(product));
192 }
193 
Log< level::Info, true > LogVerbatim
std::unique_ptr< ROOT::Math::Rotation3D > fromCMStoTB
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
static const int NCRYSTAL
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &)=0
Use this engine in event methods.
static const double deltaEta
Definition: CaloConstants.h:8
void produce(edm::Event &event, const edm::EventSetup &eventSetup) override
Produce digis out of raw data.
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
std::unique_ptr< EcalTBCrystalMap > theTestMap
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
HLT enums.
const edm::EDGetTokenT< edm::HepMCProduct > GenVtxToken
def move(src, dest)
Definition: eostools.py:511
Definition: event.py:1
#define LogDebug(id)
EcalTBMCInfoProducer(const edm::ParameterSet &ps)
Constructor.
~EcalTBMCInfoProducer() override=default
Destructor.