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