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