CMS 3D CMS Logo
Hector.cc
Go to the documentation of this file.
1 #include "FWCore/MessageLogger/interface/MessageLogger.h"
2 #include "FWCore/ParameterSet/interface/FileInPath.h"
3 
4 #include "SimTransport/HectorProducer/interface/Hector.h"
5 
6 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
7 #include "CLHEP/Units/GlobalSystemOfUnits.h"
8 #include "CLHEP/Units/GlobalPhysicalConstants.h"
9 #include "HepMC/SimpleVector.h"
10 
11 #include "TRandom3.h"
12 
13 #include "H_Parameters.h"
14 
15 #include <cmath>
16 
17 Hector::Hector(const edm::ParameterSet& param,
18  const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord>& token,
19  bool verbosity,
20  bool FP420Transport,
21  bool ZDCTransport)
22  : tok_pdt_(token), m_verbosity(verbosity), m_FP420Transport(FP420Transport), m_ZDCTransport(ZDCTransport) {
23  // Create LHC beam line
24  edm::ParameterSet hector_par = param.getParameter<edm::ParameterSet>("Hector");
25 
26  // User definitons
27  lengthfp420 = hector_par.getParameter<double>("BeamLineLengthFP420");
28  m_rpp420_f = (float)hector_par.getParameter<double>("RP420f");
29  m_rpp420_b = (float)hector_par.getParameter<double>("RP420b");
30  lengthzdc = hector_par.getParameter<double>("BeamLineLengthZDC");
31  lengthd1 = hector_par.getParameter<double>("BeamLineLengthD1");
32  beam1filename = hector_par.getParameter<string>("Beam1");
33  beam2filename = hector_par.getParameter<string>("Beam2");
34  m_rppzdc = (float)lengthzdc;
35  m_rppd1 = (float)lengthd1;
36  m_smearAng = hector_par.getParameter<bool>("smearAng");
37  m_sigmaSTX = hector_par.getParameter<double>("sigmaSTX");
38  m_sigmaSTY = hector_par.getParameter<double>("sigmaSTY");
39  m_smearE = hector_par.getParameter<bool>("smearEnergy");
40  m_sig_e = hector_par.getParameter<double>("sigmaEnergy");
41  etacut = hector_par.getParameter<double>("EtaCutForHector");
42 
43  theCorrespondenceMap.clear();
44 
45  if (m_verbosity) {
46  edm::LogInfo("HectorSetup") << "===================================================================\n"
47  << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
48  << " * * \n"
49  << " * --<--<-- A fast simulator --<--<-- * \n"
50  << " * | --<--<-- of particle --<--<-- * \n"
51  << " * ----HECTOR----< * \n"
52  << " * | -->-->-- transport through-->-->-- * \n"
53  << " * -->-->-- generic beamlines -->-->-- * \n"
54  << " * * \n"
55  << " * JINST 2:P09005 (2007) * \n"
56  << " * X Rouby, J de Favereau, K Piotrzkowski (CP3) * \n"
57  << " * http://www.fynu.ucl.ac.be/hector.html * \n"
58  << " * * \n"
59  << " * Center for Cosmology, Particle Physics and Phenomenology * \n"
60  << " * Universite catholique de Louvain * \n"
61  << " * Louvain-la-Neuve, Belgium * \n"
62  << " * * \n"
63  << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
64  << " Hector configuration: \n"
65  << " m_FP420Transport = " << m_FP420Transport << "\n"
66  << " m_ZDCTransport = " << m_ZDCTransport << "\n"
67  << " lengthfp420 = " << lengthfp420 << "\n"
68  << " m_rpp420_f = " << m_rpp420_f << "\n"
69  << " m_rpp420_b = " << m_rpp420_b << "\n"
70  << " lengthzdc = " << lengthzdc << "\n"
71  << " lengthd1 = " << lengthd1 << "\n"
72  << "===================================================================\n";
73  }
74  edm::FileInPath b1(beam1filename.c_str());
75  edm::FileInPath b2(beam2filename.c_str());
76 
77  // construct beam line for FP420: .
78  if (m_FP420Transport && lengthfp420 > 0.) {
79  m_beamlineFP4201 = new H_BeamLine(1, lengthfp420 + 0.1); // (direction, length)
80  m_beamlineFP4202 = new H_BeamLine(-1, lengthfp420 + 0.1); //
81  m_beamlineFP4201->fill(b1.fullPath(), 1, "IP5");
82  m_beamlineFP4202->fill(b2.fullPath(), -1, "IP5");
83  m_beamlineFP4201->offsetElements(120, -0.097);
84  m_beamlineFP4202->offsetElements(120, +0.097);
85  m_beamlineFP4201->calcMatrix();
86  m_beamlineFP4202->calcMatrix();
87  } else {
88  if (m_verbosity)
89  edm::LogVerbatim("HectorSetup") << "Hector: WARNING: lengthfp420= " << lengthfp420;
90  }
91 
92  if (m_ZDCTransport && lengthzdc > 0. && lengthd1 > 0.) {
93  // construct beam line for ZDC:
94  m_beamlineZDC1 = new H_BeamLine(1, lengthzdc + 0.1); // (direction, length)
95  m_beamlineZDC2 = new H_BeamLine(-1, lengthzdc + 0.1); //
96  m_beamlineZDC1->fill(b1.fullPath(), 1, "IP5");
97  m_beamlineZDC2->fill(b2.fullPath(), -1, "IP5");
98  m_beamlineZDC1->offsetElements(120, -0.097);
99  m_beamlineZDC2->offsetElements(120, +0.097);
100  m_beamlineZDC1->calcMatrix();
101  m_beamlineZDC2->calcMatrix();
102 
103  // construct beam line for D1:
104  m_beamlineD11 = new H_BeamLine(1, lengthd1 + 0.1); // (direction, length)
105  m_beamlineD12 = new H_BeamLine(-1, lengthd1 + 0.1); //
106  m_beamlineD11->fill(b1.fullPath(), 1, "IP5");
107  m_beamlineD12->fill(b2.fullPath(), -1, "IP5");
108  m_beamlineD11->offsetElements(120, -0.097);
109  m_beamlineD12->offsetElements(120, +0.097);
110  m_beamlineD11->calcMatrix();
111  m_beamlineD12->calcMatrix();
112  } else {
113  if (m_verbosity)
114  edm::LogVerbatim("HectorSetup") << "Hector: WARNING: lengthzdc= " << lengthzdc << "lengthd1= " << lengthd1;
115  }
116 }
117 
118 Hector::~Hector() {
119  for (std::map<unsigned int, H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
120  delete (*it).second;
121  }
122 
123  delete m_beamlineFP4201;
124  delete m_beamlineFP4202;
125  delete m_beamlineZDC1;
126  delete m_beamlineZDC2;
127  delete m_beamlineD11;
128  delete m_beamlineD12;
129 }
130 
131 void Hector::clearApertureFlags() {
132  m_isStoppedfp420.clear();
133  m_isStoppedzdc.clear();
134  m_isStoppedd1.clear();
135 }
136 
137 void Hector::clear() {
138  for (std::map<unsigned int, H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
139  delete (*it).second;
140  };
141  m_beamPart.clear();
142  m_direct.clear();
143  m_eta.clear();
144  m_pdg.clear();
145  m_pz.clear();
146  m_isCharged.clear();
147 }
148 
149 void Hector::add(const HepMC::GenEvent* evt, const edm::EventSetup& iSetup) {
150  H_BeamParticle* h_p;
151  unsigned int line;
152 
153  for (HepMC::GenEvent::particle_const_iterator eventParticle = evt->particles_begin();
154  eventParticle != evt->particles_end();
155  ++eventParticle) {
156  if ((*eventParticle)->status() == 1) {
157  if (abs((*eventParticle)->momentum().eta()) > etacut) {
158  line = (*eventParticle)->barcode();
159  if (m_beamPart.find(line) == m_beamPart.end()) {
160  double charge = 1.;
161  m_isCharged[line] = false; // neutrals
162  HepMC::GenParticle* g = (*eventParticle);
163  pdt = iSetup.getHandle(tok_pdt_);
164  const ParticleData* part = pdt->particle(g->pdg_id());
165  if (part) {
166  charge = part->charge();
167  }
168  if (charge != 0)
169  m_isCharged[line] = true; //charged
170  double mass = (*eventParticle)->generatedMass();
171 
172  h_p = new H_BeamParticle(mass, charge);
173 
174  double px, py, pz;
175  px = (*eventParticle)->momentum().px();
176  py = (*eventParticle)->momentum().py();
177  pz = (*eventParticle)->momentum().pz();
178 
179  h_p->set4Momentum(px, py, pz, (*eventParticle)->momentum().e());
180 
181  // from mm to um
182  double XforPosition = (*eventParticle)->production_vertex()->position().x() / micrometer; //um
183  double YforPosition = (*eventParticle)->production_vertex()->position().y() / micrometer; //um
184  double ZforPosition = (*eventParticle)->production_vertex()->position().z() / meter; //m
185  // crossing angle (beam tilt) is not known a priory; keep now 0.0 but, in principle, can be entered as parameters
186  double TXforPosition = 0.0, TYforPosition = 0.0; //urad
187 
188  // Clears H_BeamParticle::positions and sets the initial one
189  h_p->setPosition(XforPosition, YforPosition, TXforPosition, TYforPosition, ZforPosition);
190 
191  m_beamPart[line] = h_p;
192  m_direct[line] = 0;
193  m_direct[line] = (pz > 0) ? 1 : -1;
194 
195  m_eta[line] = (*eventParticle)->momentum().eta();
196  m_pdg[line] = (*eventParticle)->pdg_id();
197  m_pz[line] = (*eventParticle)->momentum().pz();
198 
199  if (m_verbosity) {
200  edm::LogVerbatim("HectorEventProcessing")
201  << "Hector:add: barcode = " << line << " status = " << g->status() << " PDG Id = " << g->pdg_id()
202  << " mass = " << mass << " pz = " << pz << " charge = " << charge
203  << " m_isCharged[line] = " << m_isCharged[line];
204  }
205  } // if find line
206  } // if eta > 8.2
207  } // if status
208  } // for loop
209 }
210 
211 void Hector::filterFP420(TRandom3* rootEngine) {
212  unsigned int line;
213  H_BeamParticle* part;
214  std::map<unsigned int, H_BeamParticle*>::iterator it;
215 
216  bool is_stop;
217  int direction;
218 
219  float x1_420;
220  float y1_420;
221 
222  if (!m_beamPart.empty() && lengthfp420 > 0.) {
223  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
224  line = (*it).first;
225  part = (*it).second;
226 
227  if (m_verbosity) {
228  edm::LogVerbatim("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line;
229  }
230  if ((*m_isCharged.find(line)).second) {
231  direction = (*m_direct.find(line)).second;
232  if (m_smearAng) {
233  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
234  if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
235  part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
236  } else {
237  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
238  part->smearAng(STX, STY, rootEngine);
239  }
240  }
241  if (m_smearE) {
242  if (m_sig_e) {
243  part->smearE(m_sig_e, rootEngine);
244  } else {
245  part->smearE(SBE, rootEngine); // in GeV, default is SBE=0.79
246  }
247  }
248  if (direction == 1 && m_beamlineFP4201 != nullptr) {
249  part->computePath(m_beamlineFP4201);
250  is_stop = part->stopped(m_beamlineFP4201);
251  if (m_verbosity)
252  edm::LogVerbatim("HectorEventProcessing")
253  << "Hector:filterFP420: barcode = " << line << " positive is_stop= " << is_stop;
254  } else if (direction == -1 && m_beamlineFP4202 != nullptr) {
255  part->computePath(m_beamlineFP4202);
256  is_stop = part->stopped(m_beamlineFP4202);
257  if (m_verbosity)
258  edm::LogVerbatim("HectorEventProcessing")
259  << "Hector:filterFP420: barcode = " << line << " negative is_stop= " << is_stop;
260  } else {
261  is_stop = true;
262  if (m_verbosity)
263  edm::LogVerbatim("HectorEventProcessing")
264  << "Hector:filterFP420: barcode = " << line << " 0 is_stop= " << is_stop;
265  }
266 
267  //propagating
268  m_isStoppedfp420[line] = is_stop;
269  if (m_verbosity)
270  edm::LogVerbatim("HectorEventProcessing")
271  << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
272 
273  if (!is_stop) {
274  if (direction == 1)
275  part->propagate(m_rpp420_f);
276  if (direction == -1)
277  part->propagate(m_rpp420_b);
278  x1_420 = part->getX();
279  y1_420 = part->getY();
280  if (m_verbosity)
281  edm::LogVerbatim("HectorEventProcessing")
282  << "Hector:filterFP420: barcode = " << line << " x= " << x1_420 << " y= " << y1_420;
283 
284  m_xAtTrPoint[line] = x1_420;
285  m_yAtTrPoint[line] = y1_420;
286  m_TxAtTrPoint[line] = part->getTX();
287  m_TyAtTrPoint[line] = part->getTY();
288  m_eAtTrPoint[line] = part->getE();
289  }
290  } // if isCharged
291  else {
292  m_isStoppedfp420[line] = true; // imply that neutral particles stopped to reach 420m
293  if (m_verbosity)
294  edm::LogVerbatim("HectorEventProcessing")
295  << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
296  }
297 
298  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
299  } // if ( m_beamPart.size() )
300 }
301 
302 void Hector::filterZDC(TRandom3* rootEngine) {
303  unsigned int line;
304  H_BeamParticle* part;
305  std::map<unsigned int, H_BeamParticle*>::iterator it;
306 
307  bool is_stop_zdc = false;
308  int direction;
309 
310  if (!m_beamPart.empty() && lengthzdc > 0.) {
311  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
312  line = (*it).first;
313  part = (*it).second;
314  if (m_verbosity) {
315  edm::LogVerbatim("HectorEventProcessing")
316  << "Hector:filterZDC: barcode = " << line << " charge = " << (*m_isCharged.find(line)).second;
317  if (m_FP420Transport)
318  edm::LogVerbatim("HectorEventProcessing") << " isStoppedFP420 =" << (*m_isStoppedfp420.find(line)).second;
319  }
320  // if ( ((*m_isStoppedfp420.find(line)).second) && ((*m_isCharged.find(line)).second) ) {
321  if (((*m_isCharged.find(line)).second)) {
322  if (m_verbosity)
323  edm::LogVerbatim("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " propagated ";
324 
325  direction = (*m_direct.find(line)).second;
326  if (m_verbosity)
327  edm::LogVerbatim("HectorEventProcessing")
328  << "Hector:filterZDC: barcode = " << line << " direction = " << direction;
329  if (m_smearAng) {
330  if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
331  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
332  part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
333  } else {
334  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
335  part->smearAng(STX, STY, rootEngine);
336  }
337  }
338  if (m_smearE) {
339  if (m_sig_e) {
340  part->smearE(m_sig_e, rootEngine);
341  } else {
342  part->smearE(SBE, rootEngine); // in GeV, default is SBE=0.79
343  }
344  }
345  if (direction == 1 && m_beamlineZDC1 != nullptr) {
346  part->computePath(m_beamlineZDC1);
347  is_stop_zdc = part->stopped(m_beamlineZDC1);
348  m_isStoppedzdc[line] = is_stop_zdc;
349  if (m_verbosity)
350  edm::LogVerbatim("HectorEventProcessing")
351  << "Hector:filterZDC: barcode " << line << " positive is_stop_zdc= " << is_stop_zdc;
352  } else if (direction == -1 && m_beamlineZDC2 != nullptr) {
353  part->computePath(m_beamlineZDC2);
354  is_stop_zdc = part->stopped(m_beamlineZDC2);
355  m_isStoppedzdc[line] = is_stop_zdc;
356  if (m_verbosity)
357  edm::LogVerbatim("HectorEventProcessing")
358  << "Hector:filterZDC: barcode " << line << " negative is_stop_zdc= " << is_stop_zdc;
359  } else {
360  m_isStoppedzdc[line] = true;
361  if (m_verbosity)
362  edm::LogVerbatim("HectorEventProcessing")
363  << "Hector:filterZDC: barcode " << line << " 0 is_stop_zdc= " << is_stop_zdc;
364  }
365  }
366  // if stopfp420 charged particles
367  else {
368  m_isStoppedzdc[line] = true; // neutrals particles considered as stopped in propagating to ZDC
369  if (m_verbosity)
370  edm::LogVerbatim("HectorEventProcessing")
371  << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
372  }
373 
374  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
375  } // if ( m_beamPart.size() )
376 }
377 
378 void Hector::filterD1(TRandom3* rootEngine) {
379  unsigned int line;
380  H_BeamParticle* part;
381  std::map<unsigned int, H_BeamParticle*>::iterator it;
382 
383  bool is_stop_d1;
384  int direction;
385 
386  float x1_d1;
387  float y1_d1;
388 
389  if (!m_beamPart.empty() && lengthd1 > 0.) {
390  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
391  line = (*it).first;
392  part = (*it).second;
393  if (m_verbosity)
394  edm::LogVerbatim("HectorEventProcessing")
395  << "Hector:filterD1: barcode = " << line << " isStoppedZDC =" << (*m_isStoppedzdc.find(line)).second;
396  if (((*m_isStoppedzdc.find(line)).second) || !((*m_isCharged.find(line)).second)) {
397  if (m_verbosity)
398  edm::LogVerbatim("HectorEventProcessing") << "Hector:filterD1: barcode = " << line << " propagated ";
399 
400  direction = (*m_direct.find(line)).second;
401  if (m_verbosity)
402  edm::LogVerbatim("HectorEventProcessing") << "Hector:filterD1:direction=" << direction;
403  if (m_smearAng) {
404  if (m_sigmaSTX > 0. && m_sigmaSTY > 0.) {
405  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
406  part->smearAng(m_sigmaSTX, m_sigmaSTY, rootEngine);
407  } else {
408  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
409  part->smearAng(STX, STY, rootEngine);
410  }
411  }
412  if (m_smearE) {
413  if (m_sig_e) {
414  part->smearE(m_sig_e, rootEngine);
415  } else {
416  part->smearE(SBE, rootEngine); // in GeV, default is SBE=0.79
417  }
418  }
419  if (direction == 1 && m_beamlineD11 != nullptr) {
420  part->computePath(m_beamlineD11);
421  is_stop_d1 = part->stopped(m_beamlineD11);
422  m_isStoppedd1[line] = is_stop_d1;
423  if (m_verbosity)
424  edm::LogVerbatim("HectorEventProcessing")
425  << "Hector:filterD1 barcode " << line << " positive is_stop_d1 = " << is_stop_d1;
426  } else if (direction == -1 && m_beamlineD12 != nullptr) {
427  part->computePath(m_beamlineD12);
428  is_stop_d1 = part->stopped(m_beamlineD12);
429  m_isStoppedd1[line] = is_stop_d1;
430  if (m_verbosity)
431  edm::LogVerbatim("HectorEventProcessing")
432  << "Hector:filterD1 barcode " << line << " negative is_stop_d1 = " << is_stop_d1;
433  } else {
434  is_stop_d1 = true;
435  m_isStoppedd1[line] = is_stop_d1;
436  if (m_verbosity)
437  edm::LogVerbatim("HectorEventProcessing")
438  << "Hector:filterD1 barcode " << line << " 0 is_stop_d1 = " << is_stop_d1;
439  }
440  //propagating
441  if (!is_stop_d1) {
442  part->propagate(lengthd1);
443  x1_d1 = part->getX();
444  y1_d1 = part->getY();
445  m_xAtTrPoint[line] = x1_d1;
446  m_yAtTrPoint[line] = y1_d1;
447  m_TxAtTrPoint[line] = part->getTX();
448  m_TyAtTrPoint[line] = part->getTY();
449  m_eAtTrPoint[line] = part->getE();
450  }
451  } // if stopzdc
452  else {
453  m_isStoppedd1[line] = false; // not stopped in propagating to ZDC and therefore in propagation to D1 too.
454  if (m_verbosity)
455  edm::LogVerbatim("HectorEventProcessing")
456  << "Hector:filterD1: barcode = " << line << " isStopped=" << (*m_isStoppedd1.find(line)).second;
457  }
458 
459  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
460  } // if ( m_beamPart.size() )
461 }
462 
463 int Hector::getDirect(unsigned int part_n) const {
464  std::map<unsigned int, int>::const_iterator it = m_direct.find(part_n);
465  if (it != m_direct.end()) {
466  return (*it).second;
467  }
468  return 0;
469 }
470 
471 void Hector::print() const {
472  for (std::map<unsigned int, H_BeamParticle*>::const_iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
473  (*it).second->printProperties();
474  };
475 }
476 
477 HepMC::GenEvent* Hector::addPartToHepMC(HepMC::GenEvent* evt) {
478  theCorrespondenceMap.clear();
479 
480  unsigned int line;
481 
482  HepMC::GenParticle* gpart;
483  long double tx, ty, theta, fi, energy, time = 0;
484  std::map<unsigned int, H_BeamParticle*>::iterator it;
485 
486  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it) {
487  line = (*it).first;
488  if (!m_FP420Transport)
489  m_isStoppedfp420[line] = true;
490  if (!m_ZDCTransport) {
491  m_isStoppedzdc[line] = false;
492  m_isStoppedd1[line] = true;
493  }
494  if (m_verbosity) {
495  edm::LogVerbatim("HectorEventProcessing")
496  << "Hector:addPartToHepMC: barcode = " << line << "\n"
497  << "Hector:addPartToHepMC: isStoppedfp420=" << (*m_isStoppedfp420.find(line)).second << "\n"
498  << "Hector:addPartToHepMC: isStoppedzdc=" << (*m_isStoppedzdc.find(line)).second << "\n"
499  << "Hector:addPartToHepMC: isStoppedd1=" << (*m_isStoppedd1.find(line)).second;
500  }
501  if (!((*m_isStoppedfp420.find(line)).second) ||
502  (!((*m_isStoppedd1.find(line)).second) && ((*m_isStoppedzdc.find(line)).second))) {
503  gpart = evt->barcode_to_particle(line);
504  if (gpart) {
505  tx = (*m_TxAtTrPoint.find(line)).second / 1000000.;
506  ty = (*m_TyAtTrPoint.find(line)).second / 1000000.;
507  theta = sqrt((tx * tx) + (ty * ty));
508  double ddd = 0.;
509  if (!((*m_isStoppedfp420.find(line)).second)) {
510  if ((*m_direct.find(line)).second > 0) {
511  ddd = m_rpp420_f;
512  } else if ((*m_direct.find(line)).second < 0) {
513  ddd = m_rpp420_b;
514  theta = pi - theta;
515  }
516  } else {
517  ddd = lengthd1;
518  if ((*m_direct.find(line)).second < 0)
519  theta = pi - theta;
520  }
521 
522  fi = std::atan2(tx, ty); // tx, ty never == 0?
523  energy = (*m_eAtTrPoint.find(line)).second;
524 
525  time = (ddd * meter - gpart->production_vertex()->position().z() * mm); // mm
526 
527  if (ddd != 0.) {
528  if (m_verbosity) {
529  edm::LogVerbatim("HectorEventProcessing")
530  << "Hector:: x= " << (*(m_xAtTrPoint.find(line))).second * 0.001 << "\n"
531  << "Hector:: y= " << (*(m_yAtTrPoint.find(line))).second * 0.001 << "\n"
532  << "Hector:: z= " << ddd * (*(m_direct.find(line))).second * 1000. << "\n"
533  << "Hector:: t= " << time;
534  }
535 
536  HepMC::GenVertex* vert = new HepMC::GenVertex(HepMC::FourVector((*(m_xAtTrPoint.find(line))).second * 0.001,
537  (*(m_yAtTrPoint.find(line))).second * 0.001,
538  ddd * (*(m_direct.find(line))).second * 1000.,
539  time + .001 * time));
540 
541  gpart->set_status(2);
542  vert->add_particle_in(gpart);
543  vert->add_particle_out(new HepMC::GenParticle(HepMC::FourVector(energy * std::sin(theta) * std::sin(fi),
544  energy * std::sin(theta) * std::cos(fi),
545  energy * std::cos(theta),
546  energy),
547  gpart->pdg_id(),
548  1,
549  gpart->flow()));
550  evt->add_vertex(vert);
551 
552  int ingoing = (*vert->particles_in_const_begin())->barcode();
553  int outgoing = (*vert->particles_out_const_begin())->barcode();
554  LHCTransportLink theLink(ingoing, outgoing);
555  if (m_verbosity)
556  edm::LogVerbatim("HectorEventProcessing") << "Hector:addPartToHepMC: LHCTransportLink " << theLink;
557  theCorrespondenceMap.push_back(theLink);
558 
559  if (m_verbosity)
560  edm::LogVerbatim("HectorEventProcessing")
561  << "Hector::TRANSPORTED pz= " << gpart->momentum().pz() << " eta= " << gpart->momentum().eta()
562  << " status= " << gpart->status();
563 
564  } // ddd
565  } // if gpart
566  } // if !isStopped
567 
568  else {
569  gpart = evt->barcode_to_particle(line);
570  if (gpart) {
571  // HepMC::GenVertex * vert= new HepMC::GenVertex();
572  gpart->set_status(2);
573  // vert->add_particle_in( gpart );
574  // vert->add_particle_out( gpart );
575  // evt->add_vertex( vert );
576  if (m_verbosity)
577  edm::LogVerbatim("HectorEventProcessing")
578  << "Hector::NON-transp. pz= " << gpart->momentum().pz() << " eta= " << gpart->momentum().eta()
579  << " status= " << gpart->status();
580  }
581  }
582 
583  } //for
584 
585  return evt;
586 }
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
Hector::m_direct
std::map< unsigned int, int > m_direct
Definition: Hector.h:103
b1
weight_default_t b1[25]
Definition: b1.h:9
Hector::m_isCharged
std::map< unsigned int, bool > m_isCharged
Definition: Hector.h:116
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
HectorTransport_cfi.ZDCTransport
ZDCTransport
HepMC source to be processed.
Definition: HectorTransport_cfi.py:6
Hector::m_rpp420_f
float m_rpp420_f
Definition: Hector.h:83
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Hector::m_ZDCTransport
bool m_ZDCTransport
Definition: Hector.h:123
Hector::m_eta
std::map< unsigned int, double > m_eta
Definition: Hector.h:113
Hector::m_rpp420_b
float m_rpp420_b
Definition: Hector.h:84
Hector::m_isStoppedzdc
std::map< unsigned int, bool > m_isStoppedzdc
Definition: Hector.h:105
Hector::m_verbosity
bool m_verbosity
Definition: Hector.h:121
Hector::m_sig_e
double m_sig_e
Definition: Hector.h:78
Hector::m_beamlineZDC1
H_BeamLine * m_beamlineZDC1
Definition: Hector.h:93
Hector::pdt
edm::ESHandle< ParticleDataTable > pdt
Definition: Hector.h:88
Hector::m_beamlineD11
H_BeamLine * m_beamlineD11
Definition: Hector.h:95
Hector::m_isStoppedd1
std::map< unsigned int, bool > m_isStoppedd1
Definition: Hector.h:106
Hector::m_isStoppedfp420
std::map< unsigned int, bool > m_isStoppedfp420
Definition: Hector.h:104
Hector::clear
void clear()
Definition: Hector.cc:137
Hector::m_sigmaSTX
double m_sigmaSTX
Definition: Hector.h:80
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
hcalRecHitTable_cff.energy
energy
Definition: hcalRecHitTable_cff.py:13
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
pi
const Double_t pi
Definition: trackSplitPlot.h:36
HectorTransport_cfi.FP420Transport
FP420Transport
main flag to set transport for ZDC
Definition: HectorTransport_cfi.py:7
Hector::m_eAtTrPoint
std::map< unsigned int, double > m_eAtTrPoint
Definition: Hector.h:111
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:222
Hector::beam2filename
std::string beam2filename
Definition: Hector.h:119
Hector::m_rppd1
float m_rppd1
Definition: Hector.h:86
Hector::m_beamlineZDC2
H_BeamLine * m_beamlineZDC2
Definition: Hector.h:94
nano_mu_digi_cff.float
float
Definition: nano_mu_digi_cff.py:13
Hector::filterD1
void filterD1(TRandom3 *)
Definition: Hector.cc:378
Hector::m_pz
std::map< unsigned int, double > m_pz
Definition: Hector.h:115
Hector::lengthfp420
double lengthfp420
Definition: Hector.h:72
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
mps_splice.line
line
Definition: mps_splice.py:76
Hector::m_TyAtTrPoint
std::map< unsigned int, double > m_TyAtTrPoint
Definition: Hector.h:110
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Hector::filterFP420
void filterFP420(TRandom3 *)
Definition: Hector.cc:211
Hector::theCorrespondenceMap
std::vector< LHCTransportLink > theCorrespondenceMap
Definition: Hector.h:125
Hector::Hector
Hector(const edm::ParameterSet &ps, const edm::ESGetToken< HepPDT::ParticleDataTable, PDTRecord > &, bool verbosity, bool FP420Transport, bool ZDCTransport)
Definition: Hector.cc:17
b2
weight_default_t b2[10]
Definition: b2.h:9
Hector::m_sigmaSTY
double m_sigmaSTY
Definition: Hector.h:81
Hector::print
void print() const
Definition: Hector.cc:471
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
verbosity
const int verbosity
Definition: DQMFileLoader.cc:24
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131
Hector::~Hector
virtual ~Hector()
Definition: Hector.cc:118
Hector::getDirect
int getDirect(unsigned int part_n) const
Definition: Hector.cc:463
Hector::filterZDC
void filterZDC(TRandom3 *)
Definition: Hector.cc:302
Hector::m_beamlineFP4201
H_BeamLine * m_beamlineFP4201
Definition: Hector.h:91
part
part
Definition: HCALResponse.h:20
Hector::m_smearE
bool m_smearE
Definition: Hector.h:79
Hector::etacut
double etacut
Definition: Hector.h:76
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
Hector::tok_pdt_
const edm::ESGetToken< HepPDT::ParticleDataTable, PDTRecord > tok_pdt_
Definition: Hector.h:69
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
Hector::beam1filename
std::string beam1filename
Definition: Hector.h:118
Hector::add
void add(const HepMC::GenEvent *ev, const edm::EventSetup &es)
Definition: Hector.cc:149
Hector::m_smearAng
bool m_smearAng
Definition: Hector.h:77
Hector::m_beamlineD12
H_BeamLine * m_beamlineD12
Definition: Hector.h:96
Hector::m_beamlineFP4202
H_BeamLine * m_beamlineFP4202
Definition: Hector.h:92
Hector::m_rppzdc
float m_rppzdc
Definition: Hector.h:85
Hector::addPartToHepMC
HepMC::GenEvent * addPartToHepMC(HepMC::GenEvent *event)
Definition: Hector.cc:477
Hector::lengthzdc
double lengthzdc
Definition: Hector.h:73
Hector::m_yAtTrPoint
std::map< unsigned int, double > m_yAtTrPoint
Definition: Hector.h:108
Hector::m_pdg
std::map< unsigned int, int > m_pdg
Definition: Hector.h:114
Hector::m_TxAtTrPoint
std::map< unsigned int, double > m_TxAtTrPoint
Definition: Hector.h:109
Hector::m_xAtTrPoint
std::map< unsigned int, double > m_xAtTrPoint
Definition: Hector.h:107
Hector::lengthd1
double lengthd1
Definition: Hector.h:74
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
Hector::m_beamPart
std::map< unsigned int, H_BeamParticle * > m_beamPart
Definition: Hector.h:102
Hector::m_FP420Transport
bool m_FP420Transport
Definition: Hector.h:122
Hector::clearApertureFlags
void clearApertureFlags()
Definition: Hector.cc:131
Generated for CMSSW Reference Manual by doxygen 1.8.14