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